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I. INTRODUCTION AND OVERVIEW 


System identification technology has been used 
successfully for many vehicles. Because of their 
large number of degrees of freedom and complex aero- 
dynamic interactions, the rotorcraft have always 
presented a special challenge to system identification 
methods. A completely integrated methodology has been 
developed under this NASA contract to solve this 
difficult problem. This methodology has also been 
translated into a user oriented series of computer 
programs. This volume provides basic guidelines for 
efficient and effective use of one of these computer 
programs. 

Figure 1 shows a schematic flowchart of the overall 
data processing technique for rotorcraft. The first 
step in this procedure is state estimation and instrument 
calibration. This is implemented by the computer program 
DEKFIS (for Discrete Extended Kalman Filter and Smoother) 
which implements an extended Kalman filter/sraoother using 
the Friedl and -Duffy formulation. Instrument biases and 
scale factors are estimated at this stage together with 
any state which is not measured directly. The second 
step involves estimation of the mathematical model of 
various forces, moments and interchanges. This is 
implemented in OSR (Optimal Subset Regression) computer 
program which uses a regression technique. Accurate 
estimates of parameters are obtained in the final step. 

One of two computer programs is used for this purpose. 
SCIDNT implements the maximum likelihood method for linear 
systems and NLSC IDNT extends the method to nonlinear 
rotorcraft models. 
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Figure 1 Integrated Rotorcraft System 
Identification Procedure 

Accuracy of parameter estimates may be improved 
by using flight test inputs based on the input design 
program, INDES. 

This user's manual describes the NLSCINDT 
computer program. The details of the theory and the 
particular implementation used are given in the final 
report . * 


* Hall, W.E., Gupta, N.K., Hansen, R. and Bohn, J., 
"State Estimation and Parameter Identification 
for Rotorcraft," final report on contract NASI-14549, 
May 1973. 
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This user's guide describes the structure of the NLSCIDNT 
program, the program's input and output, and illustrates them 
with example runs. It also gives guidelines for the program's 
effective use and information on selected topics which the 
reader will find helpful in using the program. 


The contract research effort which led to the results in this report was 
financially supported by the Structures Laboratory, USARTL (AVRADCOM), 
NASA Langley Research Center and NASA Ames Research Center. 
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CHAPTER II 
BACKGROUND 


The nonlinear, maximum likelihood, parameter identification 
computer program (NLSCIDNT) described in this user's guide was 
'^^■bten by Systems Control, Inc. (Vt) to evaluate rotorcraft 
stability and control coefficients from flight test data. 

The program has the following features : 

• The optimal estimates of the parameters (stability and 
control coefficients) are determined (identified) by 
minimizing the negative log likelihood cost function. 
These maximum likelihood estimates are asymptotically 
unbiased, consistent, and efficient (see Appendix A) 
[ 1 * 2 ]# 

• The minimization technique is the Levenberg - Marquardt 
method, which behaves like the steepest descent method 
when it is far from the minimum and behaves like the 
modified Newton-Raphson method when it is nearer the 
minimum. Hence, the technique becomes quadrat ically 
convergent as it nears the minimum, and at the same 
time it avoids the divergence problems often associ- 
ated with quadratic techniques when far from the 
minimum. 

• Twenty one states and 40 measurement variables are 
modeled, and any subset may be selected. States which 
are not integrated may be fixed at an input value, or 
time history data may be substituted for the state in 
the equations of motion. 

• Any aerodynamic coefficient may be expressed as a non- 
linear polynomial function of selected "expansion 
variables." This feature gives the user great flex- 
ibility to model nonlinear aerodynamics through the 
program's input without any changes in the program's 
code. Up to five expansion variables may be selected 
from a list of 17. 

■ ® The states and state sensitivities (partial derivatives 
of the states with respect to the identified parameters) 
are propagated by an efficient implementation of a 
variable-order, variable-stepsize , Adams -Bashforth pre- 
dictor and Adams -Moulton corrector algorithm due to 
Shampine and Gordon [3]. 


* 
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CHAPTER III 
PROGRAM DESCRIPTION 


3.1 PROGRAM STRUCTURE 

The overall logic of the NLSCIDNT program is presented in 
Fig. 3.1. It is intended to give the user a general overview 
of the program and, therefore, it shows only the major routines 
with a brief description of their purpose. 

The subroutine calling structure is given in Fig. 3.2. The 
user may find this figure useful in further understanding the 
program's flow and in constructing an overlay structure if one 
is needed. The program was developed on a CDC 7600 and re- 
quired 303 , 771 (octal) , 100,345 (decimal) storage locations 
(words) . 

The functions of the various subprograms are sketeched 
below: 

DRIVER is the main routine. It performs no calculations 
itself, but rather calls three routines -- INPUT, SMAIN, and 0UTPUT 
which accomplish the computational and input/output tasks. 
DRIVER'S most important function is the dimensioning of the major 
arrays and setting up most of the labeled commons in the program. 

BL0CKD is a block data routine. Its major function is to 
initialize arrays which contain labels for the measurement, con- 
trol, and expansion variables for input/output purposes. 

INPUT handles all of the program's input, checks for errors 
in the inputs, and initializes various flags within the program. 
It also prints this information as a record for the user. 

FLAG IT is used by INPUT to set pointer arrays to the state, 
measurement, and expansion variables which are to be used in the 
run . 
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Figure 3.2 Subroutine Calling Structure 


1 Entry point in STATE 
* Entry point in MEAS 




INREAD is called by INPUT to read the time histories of 
the measurement and control variables . 

PRNTP is called by INPUT to print the initial parameter 
values which were read into the program. 

SMAIN performs the main computational tasks. It searches 
for the optimal parameter estimates by means of the Levenberg- 
Marquardt minimization technique. It constrains the parameter 
values within bounds set by the user, and it automatically 
restricts the search for estimates of parameters having low 
identif iability . This latter function is accomplished by adap- 
tively modifying the effective rank of the information matrix. 

CVINIT is called by SMAIN to compute the covariances of the 
measurements and controls. 

SYMW1 is called by SMAIN to compute the eigenvalues and 
eigenvectors of the (symmetric) information matrix. It is opti- 
mized for use on symmetric matrices. The eigensystem is calcu- 
lated by the QR algorithm after the matrix has been put in tri- 
diagonal form by a Householder transformation. 

UPDATE ; (1) propagates the state estimates (x) and state 

sensitivities (3x/3§) with respect to the parameter estimates and 
computes the measurement estimates (y) and their sensitivities 
(3y/3§) by calls to appropriate subroutines; (2) calculates the 
cost function and its first and second gradients for use by SMAIN; 
and (3) stores the measurement estimate time histories for use by 
0UTPUT . 

INTGR3 is called by UPDATE to oversee the integration of the 
first-order differential equations for x and 3x/3§ forward 
one sample time increment. This is accomplished by a variable- 
order, variable-stepsize, predictor-corrector algorithm. The 
algorithm chooses its own order and stepsize to most efficiently 
satisfy relative and absolute error bounds on x and 3x/3§ set 
by the user. The resulting stepsize may be mors as well as less 
than the sample time increment. 
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STEP is called by INTGR8 to perform the actual integration 
of the equations. In general, the integration steosice is not 
a rational fraction of the sample time interval. 

INTRP is called by INTGR8 to interpolate the values of x, 
3x/30 and their time derivatives at the data sampling times. 

DXDTNL is called usually by STEP but initially by INTGR8 
to obtain the time derivatives of x and also Sx/39^ where § 
is an identified parameter. 

YAHAT is called by UPDATE to evaluate y and 3y/38 at 
the data sample times given x and 3x/30 at those times. 

STATE has. five ENTRY points. A call to STATE initializes 
variables used elsewhere in the subroutine. STAT I C sets the 
state initial condition- estimates , x(0). STAICI sets the state 

sensitivity initial condition estimates, 3x(0)/33. A call to 
XD0T computes the time derivative of the state estimates,’ dx/dt. 
A call to XTHD0T computes the .time derivative of the state sensi- 
tivities with respect to a specified parameter, d(3x/3§ i )/dt . 
Therefore, the vehicle equations of motion appear in this sub- 
routine . 

MEAS has three ENTRY points. A call to MEAS evaluates the 
measurement noise covariance matrix, R. HC0MP computes the 
measurement estimates, y. HTHCMP computes the measurement sensi 
tivities with respect to a specified parameter, 3y/3§ i - There- 
fore, the measurement instrument equations appear in this routine 

C0EF evaluates the aerodynamic coefficients required by the 
vehicle equations of motion in subroutine STATE. 

DC0EF evaluates the partial derivatives of the aerodynamic 
coefficients with respect to each of the identified parameters. 

C0NTRL finds the values of the control inputs at any speci- 
fied time by linear interpolation between the nearest sample time 
points. 
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0UIPUI performs three tasks: CD it computes the combined 

optimal parameters if the a priori information matrix is not zero 
(see Section 5,3); (2) it writes the time histories of the con- 
trols, measurements, and measurement estimates to tape; and (3) it 
uses PRTPLT to produce plots of the estimated measurement time 
history fit to the actual measurement time histories and to plot 
the control input time histories. 


In addition to these major subroutines, there are ten utility 
routines of which use is made by several of the above routines. 

ADD finds the sum of two matrices. 

EQUATE sets one matrix equal to another. 

INV finds the inverse and determinant of a matrix. 

MUDT finds the product of two matrices. 

PRNT writes a matrix on the printer. 

^DTITL reads a card containing the run title. 

SU3T finds the difference of two matrices. 

ZER0 - sets a matrix equal to the zero matrix. 

ASPERR produces a "walk-back" when an error in computation 
is detected. 


LNCNT prints the run's title at the top of each page 
PRTPLT produces plots on the printer. 


3.2 EQUATIONS OF MOTION 


The maximum likelihood identification algorithm may be 
applied to any dynamic system. NLSCIDNT was coded specifically 
to model the motions of a rotorcraft. These equations are 
written in terms of the 23 state variables and eight con- 
trols listed in Table 3.1. 
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Table 3.1 

State, Control, and Measurement Variables 


STATE VARIABLES 

INDEX 

SYMBOL 

DEFINITION 

UNITS 

1 

u 

longitudinal component of velocity 

ft/sec 

2 

V 

lateral component of velocity 

ft/sec 

3 

w 

vertical component of velocity 

ft/sec 

4 

p 

body roll rate 

rad/sec 

5 

q 

body pitch rate 

rad/sec 

6 

. 

r 

body yaw rate 

rad/sec 

7 

<i> 

Euler roll angle 

rad 

8 

e 

Euler pitch angle 

rad 

9 


Euler yaw angle 

rad 

10 

K 

collective flap angular rate 

rad/sec 

11 

he 

longitudinal flap angular rate 

rad/sec 

12 


lateral flap angular rate 

rad/sec 

13 

8 0 

collective flap angle 

rad 

14 

8 lc 

longitudinal flap angle 

rad 

15 

8 ls 

lateral flap angle 

rad 

16 

• 

V 

collective lag angular rate 

rad/sec 

17 

he : 

longitudinal lag angular rate 

rad/sec 

18 

hs 

lateral lag angular rate 

rad/sec 

19 

5 9r 

rotor speed variation 

rad/sec 

20 

>o 

collective lag angle 

rad 

21 

he 

longitudinal lag angle 

rad 

22 

hs 

lateral lag angle 

rad 

23 

hr 

rotor angular position variation 

rad 
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Table 3.1 (Continued) 


CONTROL VARIABLES 



SYMBOL 

DEFINITION 

6 o 

collective pitch of blades 

CD 
►— * 
n 

lateral cyclic pitch of blades 

9 ls 

longitudinal cyclic pitch of blades 

5 TR 

pitch of tail rotor blades 

6 e 

elevator angle 

6 

a 

aileron angle 

6 

r 

rudder angle 

4 f 

flaperon angle 


UNITS 



MEASUREMENT VARIABLES 




SYMBOL 

DEFINITION 

UNITS 

a x 

longitudinal accelerometer 

ft/sec" 

a y 

lateral accelerometer 

ft/ sec 2 
2 

a z 

vertical accelerometer 

ft/sec 

P m 

roll angular accelerometer 

rad/sec 2 

P m 

pitch angular accelerometer 

rad/sec 


yaw angular accelerometer 

rad/sec 2 

p m 

roll rate gyro 

rad/sec 

Pm 

pitch rate gyro 

rad/sec 

r m 

yaw rate gyro 

rad/sec 

^m 

roll position gyro 

rad 

9 m 

pitch position gyro 

rad 


yaw position gyro 

rad 

a m 

angle-of-attack vane 

rad 

6 m 

sideslip vane 

rad 

V m - 

pitot tube 

ft/sec 













Table 3.1 (Continued) 


MEASUREMENT VARIABLES ( Coni' d) 


INDEX 

SYMBOL 

DEFINITION 

UNITS 

16 

u m 

longitudinal velocity 

ft/sec 

17 

v m 

lateral velocity 

ft/ sec 

18 

w m 

vertical velocity 

ft/ sec 

.19 


rotor longitudinal force 

lbs 

20 

% 

rotor lateral force 

lbs 

21 


rotor vertical force 

lbs 

22 

^*Rm 

rotor roll moment 

ft-lbs 

23 

M Rm 

rotor pitch moment 

ft-lbs 

24 

N Rra 

rotor yaw moment 

ft-1 bs 

25 

8 1», 

flap angle of blade 1 

rad 

26 


flap angle of blade 2 

rad 

27 

S 3m 

flap angle of blade 3 

rad 

28 

B 4m 

flap angle of blade 4 

rad 

29 

B 5m 

flap angle of blade 5 

rad 

30 

6 6m 

flap angle of blade 6 

rad 

31 

6 7m 

flap angle of blade 7 

rad 

32 

^lm 

lag angle of blade 1 

rad 

33 

? 2m 

lag angle of blade 2 

rad 

34 

^3m 

lag angle of blade 3 

rad 

35 

^4m 

lag angle of blade 4 

rad 

36 

^5m 

lag angle of blade 5 

rad 

37 

^6m 

lag angle of blade 6 

rad 

38 

c 7m 

lag angle of blade 7 

rad 

39 

(cos Vm 

cosine of rotor azimuth 

N.D. 

40 

(sin Vm 

sine of rotor azimuth 

N.D. 


L3 










The equations of motion are lengthy for the sake of com- 
pleteness. As a result, they contain terms which will frequently 
be negligible, and they contain far more aerodynamic coefficients 
than can be accurately identified with ordinary rotorcraft 
instrumentation. With properly designed control inputs and with 
highly accurate and very complete instrumentation such that all 
state variables are redundantly observable, it would be possible 
to use the complete set of equations of motion and identify all 
significant parameters. However, in the common use of the pro- 
gram, only subsets of the equations are integrated and only the 
parameters appearing in those integrated equations are identified. 

For this purpose, the user has the option of fixing any state 
at its initial value or reading in a time history for it rather 
than letting the program propagate that state by integrating 
its differential equation. A considerable savings in computation 
time is realized if states unnecessary to the problem at hand are 
not integrated. As an example of the use of these options, suppose 
lateral directional data about a trim condition is to be processed 
but significant perturbations in pitch, angle-of-attack and 
flapping occurred. Then the user may choose to have the lateral 
states v, p, r, and <p integrated; to hold u, the rotor lag and 
speed states at their initial values, and to substitute measure- 
ments of q, 9, w (perhaps derived from angle-of-attack) and 
the six flapping states for these state variables in the equations 
of motion. 

The 13 degree-of-freedom, nonlinear equations of motion are 
in body-axis coordinates with origin at the center of gravity. 

Wx = f (3.1) 

The state vector is 

x = [u v w p q r * e « s 0 B lc 3 ls S Q S lc 6 U 

«o 'lc 'is SQ R c o 'lc 'is • C 3 - 2 ) 
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The vector f is composed of the following elements: 


f u = vr • V 8 sine * 5 35,2 r4 ” { C x + T aero sinS R 


' H aero cose R l 


(3.3) 


£, f - wp - ur ♦ g cose sin* * 1 on 2 R 4 s ( C y * Y aer() } (3.4) 


m 


V s u 9 - VP + 2 cos0 cos 4. + i pfi 2 R 4 tt { C z T aerQ cos0 R 


• H aero sin0 R } 


(3.5) 


f_ = 


.2 2 , 


(I„-I_)qr + I _ (q -r ) + (I q - I r)p 
y z ^ yz n v xz^ xy v 

2 5 

+ pQ R it { C T + Q ■ . sine. - L cosS^ 
L x aero R aero R 


• Y aero Z hub + N b cose R^o^ . * 6 


6a w lc " I Ca C lc )} ^ 3 ' 6) 


f = (I -I )pr + I (r“-p Z ) + (I r - I q)p 
q z x'^xz- Vi v xy yz M ^ H 


* ° n2 RS,r fC M * M aero * H aero (x hub sin9 R z hub 005 V 


* T aero (x hub cos9 » * z w, sin8 D ) * N K [I S „ 


hub 


b L 8a w ls 


* Ha ? ls ’ ! 0 (P' cos9 r * r' sine R )]) 

f r - t I x’V Pq * I xy ( P 2 -'l 2 ’ + f I y Z P- it x2 q > r 

* ° n2 R$ * 1 C N * Y aero x hub ' L aero sinS R 


(3.7) 


^aero cos9 R * N b sine R ( ' I S a 6 l'c ‘ ? ic * I o q ' )) (3 ' 8) 


f ^ . a p + (q sin<}> + r cos<{>) tane 


(3.9) 
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i 


f 9 “ 

q cosij) - r sin4> 

(3.10) 

f * * 

(q sinp + r. cos<f>)/cos0 

(3.11) 

f S 0 ’ ' 2I Bi (p sine R • r cos8 R^ ' Vb 6 o 



+ '.pfl 2 R 5 tt M° 

(3.12) 

f 8 

B lc 

X 0^ V 8 ’ 1 - )S lc ‘ CI B S s v g + 2I 68- )B 1 s 

(3.13) 


+ 2I ga (p cos0 R + r sin0 R ) + pfi 2 R 5 it M 1c 

(3.13) 

£ *i, 

* (I 8 «s U 8 + 2I SS )8 lc : VV 1)8 ls ' 2I S a<! 



+ oil 2 R 5 tt M 1s 

(3.14) 

V 

8 o 

(3.15) 


■ 8 lc 

(3.16) 

£ hs 

■- 8 ls 

(3.17) 

11 

o 

4h 

‘ 2I ci fp sin9 R ‘ r cos9 r) * S 5 ” Q° 

(3.18) 

%c 

= -I c (v 2 -l)C lc - (I. g s * 2I c -)c ls - pa 2 R 5 »Q lc 

(3.19) 


^ 8 s v ? * 2I ?c^lc ■ I ;^' , ;' 1 - I5 ls 



2 2 Is 

+ pjr r^tt q 1s 

(3.20) 

f 6Q R 

= pfi 2 R 5 tt Q n 

(3.21) 

f, 

"o 

^0 

(3.22) 
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(3.23) 


f ‘i= ‘ 


f r s ?1 S 

g ls iS 

f «n R * S 


(3.24) 

(3.25) 


where 


C x = C x 0 * C x u OR * C x v SR * C x w Em + c x p % 


* % “ * C *r'0 + C xi TR 6 TR + C x Se a e 


+ C6 + C 6 + C 5/- 

X6 a a x«s r r ^x 6f °f 


(3.26) 


and Cy» C^, C^, and are expanded similarly; 


aero 


T o *X T* ^ * T «ls ^ + V. 


* T »leV * T »!> 7 T i 0 T * T C lc T 2 


+ T 


ls + T, 40r 


‘hi a ' ‘ sn R f! * T V° + T nc' lc 


+ T r C, + T - . 6<J/ R + T fl 0 n + T fi 0. + T fl 0- 

5 Is ls 6 ^R R e 0 0 9 lc lc 9 ls ls 


(3.27) 


and H , Y , L , M__ and Q a „_ are expanded 
aero aero aero aero* ^aero r 

similarly; 
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- M 0 ° * ♦ m;-* * < A * $ * M° §.♦ M? § 


6o B : » 


* M* TT * M 


Bis “ -Bo « 8 lc lc 

o SII R 

8 B,, ♦ Mr -,s * MT — ^ + M* C lc * M?„ 

8 ls 15 5 o a tie n Cis ls «! R n- 


. wO fl . u O ^0 . ,,0 ? lc 
+ Mq P ~ + M* “FT + M* 


+ M °o'° + M °ic' lc + M °is' ls + M ^ R 6 ^R + M e o e o 


+ M° 0 1 + M? 0. 

9 lc lc 9 ls ls 


(3 


and M^ c , M^ s , Q°, Q^ c , Q^ s , and Q— are expanded similarly. 
The elements of the matrix W which are not zero are: 


W uu ■ 1 


Vis 


_ pfl 2 R 4 TT 


m 


N b • Sin 6 R 


o n 2„4 N, S_ 

W . = - -j- -4 COS 0 D 

u 6o m 2 n 2 R 


W = 1 
W 


w . »ii¥i 

v ?lc " 


N b . s c 


W - 1 

WW 


.28) 
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w . - ££^I • N. . 4 cos e 

wg 0 m b R 


rt0 2 0 4 N K S_ 

W . = p - fl m R * -4 sin 0p 

w ? ls m 2 n 2 R 


W = I 
PP x 


W = -T 
pq A xy 


W = -I 
pr xz 


W . = -pfl 2 R 5 TT ^ • COS e 

PBis 1 R 


\2r,5. 


T • 

A Bip . 


w pg = pQ R tt . C2N b ) -g* ; sin 


,2 n 5 


Ba 


W_ = Pfl^ir N b -g=- cos 0 


P6 


lc 


>2 n 5. 


'6a 


W P6 1S " : pQ “ R-,r N b if C0S 9 R 


>2^5. 


W p . = -pQ-R-ir N b -g* sin 0 R 


W 


P^lc 


. n0 2„S ff N b £5, Z hub 
pnRir T^2 R 


W 


P^, 


•pa 2 R 5 TT — y 
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W c = pfl 2 R 5 TT (2N b ) sin 0 R 
^ o 


!• 

W_. = -pfl 2 R 5 TT-N b -I 2 - cos 0 R 


p? ls 


W = -I 
qp xy 


W = I 

qq y 


^qr ” ’ *yz 


,2 n 5_ xf S 8 Z hub 


hub 


\h - p!i R ” N b $ nr sin s r + nr cos e R 


V,. * -^ 2r5 " ^rrr 


1C 


S2 


W qB lc " N b -S 2 - 


W qB ls = - PQ2R51T N b -IT 


\l, * pn2R5lT "T (^TT cos 9 R • HT sin 9 R 

x S M 


X • 

w nr = -pfl 2 R 5 ir N, 

qc lc b S2 


w rp ' •'« 
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W a -I 
rq yz 


W = I 
rr z 


W • = -p£2 2 R 5 ir ^ sin e„ 

B ls 2 £2 2 R 


1 • 

W rg = *P fl2R5ir ( 2N b) cos 9 1 


W 


rB 


lc 


= p£2 2 R^tt -j~ sin 0 R 


w rB ls = -o 112115 * N b “I 4 sin 9 R 


w 




• = pfi 2 R 5 TT ^ cos 0 R 


£2 


.2t, 5 N b S C X hub 


T.T _ ^ rj U L, UU 

W rc, „ pn R 77 2 n 2 R 


lc 


£2 


W • = p£2 2 R^TT Nv — y cos 0 D 

r *s b ? R 


,2 d 5. 


I: 


W r ^ = -pn'R-’rr C2N b ) -&■ cos 0j 


W 


r ?ls ' •» [!2rS ’' N b ^ sin S R 


’ 1 
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W =1 
00 1 


W, , - 1 


\u * S B sin e R 


W* = -S. cos 0 — 
B 0 W B R 


W B o q = ' S B ( ‘ Z hub sin 0 R ’ x hub 


W* • = I 

B B 8 

p o p o p 


W 


S 0 6 0 = X B 8 S v 6 + 2I B6 


W* = 21 • 
8 o*s B<J> 


W* = -I. 
Bi c q Ba 


w* = I 

8 lc 8 lc 8 


W 


8 lc 8 lc = 10 8s + 2Ifi8 


W B 8 “ 2 1 s 

8 lc 8 ls 8 


W B ls P 3 * 1 Ba cos e R 


W 3 ls r = * 1 Ba sin 0 R 


£. Z 


COS 0 R ) 


I 



w* = I 

8 ls 8 ls 8 


\s 8 lc ■ - 2l8 


' 8 ls 8 ls * I(S ® s ^ * 2188 


V» 6 * 1 


W 8 fi 1 

B lc e lc 


W =1 

8 ls 8 ls 


W* = I sin 0 D 
C 0 P ?ct R 


W; ' = -I COS 0 n 

5a R 


Vo 5 


Vs ^ 


W Lc„ = *5 g S V 5 + 2I 55 


!• , = 21 • 


W* * S 
C lc v « 
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W. * -S z. . 
C lc P 5 hub 


W. _ = S_ x, . 
?lc r C hub 


W # • a I 

? lc ? lc 5 


W; 


5 lc c lc 


g S V C + Zl d 


W* = 21 

5 lc 5 ls 5 


W. ti = S, cos 0 
’Is 


C lo u ~R 


w ? ls w = s ; 5in 9 R 


W ? ls q “ S c tz hub cos e R " ’'hub 


w 


• « 


c ls c ls- 5 


W* = -21 

*ls c lc ^ 


W • = I £ _ v + 2 1 ■ • 

5lS C lS 5 S 5 « 


W* • = 1 

Vs 


W =1 
s o s o 


sin e R ) 
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where 


W 


? lc ? lc 


1 



ls ? ls 


1 


w 


Ms 


1 


g 

p 

m 

Si 

R 

I. 


xy 

X xz 

: yz 

SB 


Sc 


2 

acceleration of gravity, ft/sec 
density of the air, slug/ft l * 3 
mass of the rotorcraft, slug 
reference rotor speed, rad/sec 
rotor radius , ft 

2 2 2 

roll moment of inertia = /(y + z )dm, slug-ft 

2 7 2 

pitch moment of inertia = /(x +z )dm, slug-ft 


2 2 

yaw moment of inertia = f (x + y )dm, slug-ft 

? 

x-y product of inertia « /xy dm, slug-ft - 

2 

x-z product of inertia ■ fxz dm, slug-ft 

2 

y-z product of inertia =» fyz dm, slug-ft 

x 


N s V r 


All terms non-dimensional 


2 


l 

/ N m, dr 

o c b 


where 

N g = 1st flapping mode shape, N.D. 
N = 1st lagging mode shape, N.D. 
m^ = blade mass 

' implies division by QR. e.g. p' = p/n R 
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3.3 MEASUREMENT EQUATIONS 


Forty measurement instruments are modeled in the program, 
from which the user may select any subset. Different measure- 
ment equations are applicable depending upon the set of 
equations of motion selected. 

3.3.1 Measurements for Complete Set 

The equations which evaluate the measurement estimates when 
the complete set of equations of motion is used are: 


(1) longitudinal accelerometer 


a = u - vr + wq 


+ a 


sin 0 - (q + r jx 




+ Cpq - r )y_ a + Cpr + q)z_ CT + b 

°x L °x 

(2) lateral accelerometer. 


a = v - wp + ur - § sin $ cos 6 + (pq + r)x 

■ y 

- (r2 * pZ)y cg y * ^ - P5 z cg y * b y ■ 


(3) vertical accelerometer 


w - uq + vp 


cos 


<p cos 9 + (pr - q)x c 


+ (rq + p)y ca - (q 2 + p 2 ): r „ + b. 


eg. 


3L Z . — 


o 

© 



(4) roll angular accelerometer 


« • 


p = p + b* 
*m * p 


(5) pitch angular acceleromete] 


q * q + b* 

H m H q 


(6) yaw angular accelerometer 


r = r + b* 
m r 


(7) roll rate 


P = d + b 
*m y p 


(8) pitch rate 


q_ = q + b n 
m ^ q 


(9) yaw rate 


r_ ‘ “ r + b 
m r 


CIO) roll angle 

*m " * * b * 

(11) pitch angle 

= 6 + b a 
m o 

(12) yaw angle 


■ If * b t( 
in y 
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( 13 ) angle-of-attack vane 


a_ = ( 1 + k ) tan * 1 


m 


w - qx + dy 

Cg • ; C2 

°a 3 a 

u - rv +• q z ' 


( 14 ) angle-of-sideslip vane 


m 


(l + k 3 ) sin ' 1 


v + rx„„ - pz 

c g 3 H eg. 


where 


V. 


u-ry +qz 

°S C °SJ 


w '^ cg tpy • 
*8 S S 


v^ rzc -p: 

. G *3 g 


2 ) : 


(15) total velocity 


V m = (1 + K y ) /u 2 + v^ + w 2 + 


v 


(16) u = (1 + k ) u + b 
m u u 


(i 7 ) v m = Cl t k v )v * b y 


U 8 ) » m - (1 t k w )w t b w 


( 19 ) X R ^ * (1 k x )jrp(! 2 R 4 (T„„„ sin e D - H 


m 


aero 


aero 


(20) Y r = (1 + k v ) irpfl 2 R 4 Y + b v 

R m 7 R aero Yi 


•r - ( 1+k z )*Pft 2 R 4 (-H _ sin 0 D - T 

K m z r aero R aero 


+ o 



cos 


cos 


( 21 ) 


lO 



(22) L,^ - (1 - k LR ) 1 rp£2 2 R S CQ aer0 sin 9 R - L aer0 cos 9 R 

Y aero ^ref^ + 

(23) M Rm = C 1 *>S| R )s<=« 2R5 t M aero + H aero CX ref sin «R 

- 2 ref c °s e R ) * T aero (3E ref cos e R 

* “ref sin e R ) l * b M R 


(24) N = (l + k )soil 2 R 5 (ii ref Y aero -L r ._, sin 9, 
m K 


• <5aero cos 6 R> '* b 


N 


R 


( 25 ) 6 lm = (1+k 6 ^3 +6 ' fS lc + 6 lc ) cos(flt +&,) 

im a 1 o o ref ic ic re£ l 

.-(3, +8,_ ) sin(flt+e.)] + b ft 

is is ref 1 .1 

(26) S 2m " U*^) [ V B o ref ' (6 lc* 6 lc ref ) c °s(<lt+e 1 »40 1 ) 

-CBis^is ) sin(nt+e 1 + A0 1 ) ] + b g 
rs f 2 

(27) 63 = (l + kg )[6 0 +6 0 '( e ic +6 lc ^ cosCnt+e.+ZAej) 

m 3 ref ''ref 

'( 6 ls + 6 ls f ) sin C^t+ 8 1 + 2A0 1 ) ] + bg ^ 

(28) 0. = (l+k g ) [6 +6 - (0 1 +6, )cos(nt+0 1 +3A0. ) 

m b 4 0 °ref ic lc ref 1 1 

* ( - S ls + S ls )sin(fit+9 1 + 3A6 1 ) ] + b g 
ref d 4 
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(29) e 5 = Cl+k B HV S o / (S lc tS lc )cosC«t*e l *4Ae 1 ) 

m 5 ret ret 

-(B,_ + B, )sin(Qt+0, +4A0.) ] + bg 
iS 1 re£ 5 


(30) 6^- Cl*k S6 )t8 0 *S 0re{ -(B lc +S 1Cref )cos(8f« 1 *SA9 1 ) 

-(8, + 8, )sin(«t+0 1 +5A0 1 ) ] + bg 

ls ls ref 11 6 

(31) B n = (1+kg ) [B +6 )cos(nt+0 1 +6A0 1 ) 

'm .7 0 °ref . iC ref 1 

- (B, +8, . ) sin (ft t+0 , +6A0 .. ) ] + b R 

is ls re£ 1 i *7 


(32) ? x = Cl+k c )[ ? * 5 ' cos(:It * e 1 ) 

m 1 rei rer 

- U, + c ) sin (ftt+9 , ) ] + b 
is J "'re£ 1 ? 1 

(33) **2 ■ ( W t Ht o n„ ( -t'lc +5 lc )cos(!Jt + 6 1 »40 1 ) 

m Z ref ret 

-(?-, +C-, ) sin(ftt+0, + A9, ) ] + b 

ls ls ref 11 c 2 


(34) - ■(l*k )[c>t 0 -(C lc *4 lc )cos(Sit.9 *2AS ) 

m 3 ref ret * 

-Cc ls *c ls - )sin(^t+0 1 +2A0 1 ) ] + b ? 

ref 3 


(35) q . = (1 + k )U_ + C Q •^ir +i; lc )cos(flt+.e.+3A0 1 ) 

4 m ? 4 0 °ref lc lc ref 1 1 

■( 5 i« + Ci s Osinfflt+O^+SAB.)] + b 
ls is ref 1 1 c 4 

(36) t 5 = d*k )U 0 *« 0 „' (c 1c* c 1c )cos(nt*e *4Ae x ) 

m 5 rer rer 

-(c. +C,_ )sin(ftt+9, +4A0. ) ] + b r 
ls s ref 1 1 "5 


3 0 


(37) 


)cos (J2t+C 1 + 5A8 1 ) 


(38) 


(39) 

(40) (cos * R ) m = (1 *k c , R ) cos » R + b ci(jR 


5.4 POLYNOMIAL EXPANSION OF COEFFICIENTS 

Three hundred forty eight aerodynamic coefficients appear 
in the equations of motion in Section 3.3. These coefficients 
are listed in Table 3.2 by giving the index of the coefficient 
in the row of the appropriate primary symbol and column of the 

appropriate subscript. For example, coefficient 46 is C T and 

o 

coefficient 112 is T* . Coefficients which are usually not 

8 o 

significant have their entries enclosed in "( )", i.e., 
coefficient (3) is C 

x v 

Each of these coefficients is evaluated within the program 
as a polynomial (although in most applications the polynomial 
will consist of a single, constant term). The procedure by 
which this is accomplished will appear complex at first; but 
once it is mastered, the user may appreciate .the great 
flexibility which it gives him. 


' 6 n, ' ( 1 *\ )U ° n o r ef‘ U lc + 'lc 


'ref 


- (Ct+Ci „ )sin(fit + 9.. + 5A0.. ) ] + b 

ls ls ref 1 1 c 6 

' 7 m " (1 * k C 7 )[! o* i: o ref - (,: lc* c lc ref )cos ( !It * 9 l* 6a9 1 ) 

+ )sin(nt+9.+6A6 1 )] + b_ 

1S iS ref 11 ? 7 

(sin V. - (l*k s<>R ) sin » R * b s( , R 



Table 3.2 

Indices of the Aei’odynamic Coefficients 


PKIMAkY 

5YMHM. 


SUUSf.KIPI 


(i/) ia 
■n (ii) 
(4/) 4H 
62 (6J) 

(7/) /a 


4 

0>) 

(19) 

20 

34 

(35) 

(«) 

50 

64 

(65) 

(/9) 

80 


W K 'o 'lc|‘U 


213 

210 m 23a 

2b 1 262 263 


/lib 28/ 288 

ill IK* JIJ 

116 J 3/ Jill 

Jo I J62 J6J 


2U 21b 216 

239 240 241 

264 26S 266 


290 291 
315 316 
340 341 
36b 366 


217 I 201 
243 


93 94 

113 114 
133 134 
153 154 
1/3 174 
193 194 

219 220 
244 245 


95 96 

115 116 
135 136 
155 156 
1/b 176 
195 196 

221 222 
246 247 


267 260 269 270 271 272 


292 29J 

317 318 

342 343 

367 J61J 


244 295 
319 320 
344 345 
369 3/0 


9a 

99 

Too 

101 

102 

103 

(104) 

(105) 

l ia 

119 

120 

121 

122 

123 

( 124) 

(125) 

138 

139 

140 

141 

142 

143 

(144) 

(145) 

158 

159 

160 

161 

162 

163 

(164) 

(165) 

178 

1/9 

180 

>“> 

182 

183 

(IH4) 

(1115) 

198 

199 

200 

201 

202 

203 

(204) 

(205) 

224 

225 

226 

227 

228 

229 

(23U) 

(231) 

249 

250 

251 

252 

253 

254 

(255) 

(256) 

2/4 

275 

2/6 

277 

2/8 

279 

(200) 

(201 ) 

299 

300 

301 

302 

303 

304 

(305) 

(306) 

324 

! 325 

326 

327 

1 328 

329 

(330) 

(331) 

349 

350 

351 

352 

353 

354 

(365) 

(356) 

374 

375 

376 

377 

378 

379 

(3it0) 

(301) 


83 | (U4) | (05) | (06) | (0/) 






































Because of its complexity, the procedure is best introduced 
by example, Suppose the rolling moment coefficient is to be 
modeled as 


C L = c * C t ca-8 0 ] * C C3-6 0 ) 2 . C L 


ao 


a 6 


TR 


TR 


where S ', is a constant (not necessarily 3 at t a 0) . Let 


Pl = C L 0 - p 2 = v p - V p «‘ c i « 4 


TR 


Z 1 = a ’ Z 1 = °’ Z 2 = S ’ z 2 = 3 0 ’ z 3 = 5 TR’ 


Z, = 0 
J o 


Then, Eq. (3.69) is identical to 


P 1 + P 2^ Z 2 ' Z 2 J + P 3^ Z 2 ‘ Z 2 ^ +P 4 Z 1 Z 3 


C T = 


Pl C2 l' Z l > < 2 2' 2 2 > C z 3' 2 3 >' 

0 0 0 


* ? 2 ( 2 1 - 2 1 5° ^ 2 2 " 2 2 (2 3 - 2 3 >' 

o o 0 


+ P 3 (Z 1 ' Z l o )0 (Z 2 * Z 2 J 6 (Z 3 * Z 5 
o 0 0 

* fir'll 1 ** 2 2 " 2 2 ^ ° C2 3‘ 2 3 J 1 
0 o 0 


2 


33 



4 

2 

j-1 


I ° 

11 

( i-1 


n 




o! 


where the exponents have appropriate values. 

The form of Eq. (3,72) which emerged in the above example 
is similar to the general form used to evaluate each of the thirty 
three coefficients. Specifically, the k-th coefficient is 
computed as 


s - 



n iii 


Pj is a parameter which may be identified by the 
program, 

is an "expansion variable," 

2 . is the reference value of the i-th expansion 

^ Q 

variable, and 

j are integer exponents of the expansion variables. 

The exponents and the expansion variable reference values are in- 
puts to the program and remain constant throughout the entire run. 
The parameters are constant throughout an iteration of theoptim- 
itation search, but they may vary from one iteration to the next 
if they are identified. The initial values of the parameters are 
inputs. The expansion variables are functions of the state and 
control variables and are evaluated repeatedly by the program. 



Table 3 . 3 

Expansion Variables 


INDEX 

SYMBOL 

DESCRIPTION 

UNITS 

1 

a 

angle-of-attack 

rad 

2 

B 

angle of sideslip 

rad 

3 

P 

roll rate 

rad/sec 

4 

q 

pitch rate 

rad/sec 

5 

r 

yaw rate 

rad/sec 

6 

9 

o 

collective pitch 

rad* 

7 

9 lc 

lateral cyclic pitch 

rad* 

8 

9 ls . 

longitudinal cyclic pitch 

rad* 

9 

6 tr 

tail rotor collective pitch 

rad* 

10 

6 

e 

elevator angle 

rad* 

11 

6 

a 

aileron angle 

rad* 

12 

6 

r 

rudder angle 

rad* 

. 13 

s f 

flaperon angle 

rad* 

14 

u 

advance ratio 

N.D. 

15 

> 

u 

longitudinal velocity (normal ized) + 

N.D. 

16 

* 

V 

lateral velocity (normalized) 

N.D. 

17 

> 

w 

vertical velocity (normalized) 

N.D. 


★ 

Recommended units; actual units are determined by the input 
time histories. 


+ u" = u/nR 
v' = v/flR 
w' = w/QR 









As many as five expansion variables may be selected for a 
run. Storage and run time constraints make more than five 
impractical. All the polynomials describing the coefficients 
must be expressed in terms of the same five (or fewer) expansion 
variables. The selected variables are specified by input from 
the list of 17 possible expansion variables in Table 3.3. They 
may be chosen in any order, but once chosen, the reference 
values and the exponents must be input such that they correspond 
to the same order. Of course, the units of the reference values 
must also be consistent with the units of the corresponding 
expansion variables. 

The equations defining the expansion variables a, 3, and 


y are : 



a 

= tan ^ 

U) 

3 

= tan V 

M * w 2 ) 

u 

= OR Cu 

cose R + w sin0 R ) 

3.5 PARAMETERS 



NLSCIDNT was written with a capacity for 705 parameters for 
use in evaluating the equations of motion, the measurment equa- 
tions, and the expansion variables (see Table 3.4). This is more 
than sufficient for the present model, so not all parameters are 
used. Also, some parameters have been reserved for uses which ar 
not currently implemented but are easily added if the need arises 
Finally, some parameters may not be identified; for example, the 
aircraft mass and the air density may not be identified. 



Table 3.4 
Parameters 


INDEX 

SYMBOL 

1-400 


401 

u(0) 

402 

v(0) 

403 

w(0) 

404 

p(0) 

405 

q (0) 

406 

r(0) 

407 

*(0) 

408 

0(0) 

409 

<M0) 

410 

6 0 (°) 

411 

®1 C (°) 

412 

6 1S (°) 

413 

S.(°) 

414 

■ B lc (0) 

415 

Bi s (°) 

416 

c 0 (°) 

417 


418 

c ls (0) 

419 

fiiJp(O) 

420 

5 0 (°) 

421 


422 

Ci s (°) 

423 

«* R (0) 

426 

9 

427 

P 

428 

®R 


DESCRIPTION 


These parameters are available for defining 
the polynomial expansions of the aerodynamic 
coefficients 


Initial conditions for the state 
variables 


Gravitational acceleration (default - 
32.174 ft/sec 2 ) 

Air density 

Rotor shaft tilt in X-Z plane, (+) forward 






Table 3.4 (Continued) 




DESCRIPTION 


Rotor speed 
Rotor radius 
NC/ R 

Lock number = pacR^/I b 
2-D lift curve slope 
Mass of vehicle 

Longitudinal distance of rotor hub from 
vehicle C.G. 

Vertical distance of rotor hub from vehicle 
C.G. 


fuselage moments of inertia 

(slug-ft^) 


blade moment of inertia (slug-ft^) 
blade flap moment of inertia N.D. 
blade lag moment of inertia N.D. 


inertia terms - for definition see 
page xi. Definitions of Inertial 
Constants in the final report [10] 


Number of blades 


58 







INDEX 

SYMBI 

459 

V 0 

460 

461 

9s 

470 

R x 

471 

b X 

472 

k X 

473 

x cgx 

474 

y cgx 

475 



cgx 

479 

R y 

480 

v 

481 

k y. 

482 

x cgy 

483 

y cgy 

484 

z cgy 

488 

R z 

489 

b z 

490 

k z 

491 

x cgz 

492 

y cgz 

493 

z cgz 


Table 3.4 (Continued) 


DESCRIPTION 


J Inertia terms 

Longitudinal Accelerometer 
Noise covariance 
Bias 

Scale factor error 
X location 
Y location . 

Z location 


Lateral Accelerometer 
Noise covariance 
Bias 

Scale factor error 
X location 
Y location . 

Z location 


Vertical Accelerometer 
Noise covariance 
Bias 

Scale factor error 
X location 
Y location 
Z location 
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Table 3.4 (Continued) 


INDEX 


SYMBOL 


DESCRIPTION 


494 

495 


497 

498 

499 


503 

504 

505 


Vertical Accelerometer (Cont'd) 


Roll Angular Accelerometer 
Noise covariance 
Bias 

Scale factor error 


Pitch Angular Accelerometer 
Noise covariance 
Bias 

Scale factor error 


509 

R; 

Yaw Angular Accelerometer 
Noise covariance 

510 

r 

b* 

r 

Bias 

511 

k r 

Scale factor error 

515 

R_ 

Roll Rate Gyro 
Noise covariance 

516 

P 

b p 

Bias 

517 

k p 

Scale factor error 











Table 3.4 (Continued) 


SYMBOL 


DESCRIPTION 



Pitch Rate Gyro 


Noise covariance 


Scale factor error 



Yaw Rate Gyro 


Noise covariance 
Bias 

Scale factor error 


Roll Position Gyro 


Noise covariance 


Scale factor error 


Pitch Position Gyro 


Noise covariance 


Scale factor error 










Table 3.4 (Continued) 


INDEX 


545 

546 

547 

548 

549 


SYMBOL 


> 


DESCRIPTION 


Yaw Position Gyro 
Noise covariance 
Bias 

Scale factor error 


• 


Angle-of-Attack Vane 

551 

R 

ct 

Noise covariance 

552 . 

b 

a 

Bias 

553 

k 

0 1 

Scale factor error 

554 

x cga 

X location 

555 

y cga 

Y location 

556 

Z cga 

Z location 

557 

Sideslip Vane 

558 

r b 

Noise covariance 

559 

b B 

Bias 

560 

k 8 ' 

Scale factor error 

. 561 

x cgB 

X location 

562 

y cgB 

Y location 

563 

564 

z cgB 

Z location 
Pitot Tube 

565 

r vt 

Noise covariance 

566 

b VT 

Bias 

567 

k VT 

Scale factor error 

568 

X VT 

X location 

569 

y VT 

Y location 

570 

Z VT 

Z location 
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Table 3.4 (Continued) 


INDEX 

SYMBOL 

DESCRIPTION 



Pitot Tube (Cont'd) 

571 

Po 

Reference air density 

572 

V S 

Velocity of sound 



Longitudinal Velocity Measurement 

574 

R u 

Noise covariance 

575 

b u 

Bias 

576 

k 

Scale factor error 


u 




Lateral Velocity Measurement 

578 

V 

Noise covariance 

579 

b V 

Bias 

580 

k v 

Scale factor error 



Vertical Velocity Measurement 

582 

R w 

Noise covariance 

583 

b w 

Bias 

584 

k w 

Scale factor error 

585 



586 

x ref 

Horizontal reference distance . 

588 

*ref 

Vertical reference distance 



Rotor Longitudinal Force Measurement 

590 

r xr 

Noise covariance 

591 

Ct 

X 

Bias 

592 

k *R 

Scale factor error 










Table 3.4 (Continued) 


INDEX 

SYMBOL 

DESCRIPTION 



Rotor Lateral Force Measurement 

594 

R YR 

Noise covariance 

595 

b YR 

Bias 

596 

k YR 

Scale factor error 


' 

Rotor Vertical Force Measurement 

598 

Rz R 

Noise covariance 

599 

bz R 

Bias 

600 

k ZR 

Scale factor error 



Rotor Roll Moment Measurement 

602 

R Lr 

Noise covariance 

603 

b l-R 

Bias 

604 

k LR 

Scale factor error 



Rotor Pitch Moment Measurement 

606 

r Mr 

Noise covariance 

607 

b M R 

Bias 

608 

k M R 

Scale factor error 



Rotor Yaw Moment Measurement 

610 

r Nr 

Noise covariance 

611 

b N R 

Bias 

612 

k N R 

Scale factor error 

614 

*1 

Euler angle 

615 

A<j> 

Increment 

616 

9 1 

Euler angle 

617 

A6 

Increment 
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Table 3.4 (Continued) 


INDEX 

SYMBOL 

618 

®o ref 

619 

®lc ref 

620 

S ls ref 

622 

R B1 

623 

b 81 

624 

k Bl 

626 

R 82 

627 

b 82 

628 

k 82 

630 

R B3 

631 

b B3 

632 

k B3 

634 

R B4 

635 

b B4 

636 

k B4 

638 

R B5 

639 

b B5 

640 

k B5 


DESCRIPTION 


Blade flapping 


Blade 1 Flapping Measurement 


Noise covariance 
Bias 

Scale factor error 
Blade 2 Flapping Measurement 


Noise covariance 
Bias 

Scale factor error 
Blade 3 Flapping Measurement 


Noise covariance 
Bias 

Scale factor error 
Blade 4 Flapping Measurement 


Noise covariance 
Bias 

Scale factor error 
Blade 5 Flapping Measurement 


Noise covariance 
Bias 

Scale factor error 








Table 3.4 (Continued) 


INDEX 

SYMBOL 

642 

R B6 

643 

b 06 

644 

k B6 

645 


646 

R B7 

647 

»B7 

648 

649 

k B7 

650 

5 o ref 

651 

? lc ref 

652 

653 

^ls ref 

654' 

\l 

655 

b .l 

656 

k Sl 

657 


658 

R C2 

659 

b C2 

660 

661 


662 


663 

\z 

664 

665 

k ;3 


DESCRIPTION 

Blade 6 Flapping Measurement 
Noise covariance 
Bias 

Scale factor error 

Blade 7 Flapping Measurement 
Noise covariance 
Bias 

Scale factor error 


Lag references 

Blade 1 Lag Measurement 
Noise covariance 
Bias 

Scale factor error 

Blade 2 Lag Measurement 
Noise covariance 
Bias 

Scale factor error 

Blade 3 Lag Measurement 
Noise covariance 
Bias 

Scale factor error 








Table 3.4 (Continued) 


INDEX 

SYMBOL 

DESCRIPTION 



Blade 4 Lag Measurement 

666 

R 54 

Noise covariance 

667 

b S4 

Bias 

668 

669 

k ;4 

Scale factor error 
Blade 5 Lag Measurement 

670 

R C5 

Noise covariance 

671 

b ?5 

Bias 

672 

673 

k 55 

Scale factor error 
Blade 6 Lag Measurement 

674 

\ 6 

Noise covariance 

675 

b S6 

Bias 

676 

k S6 

Scale factor error 

677 

Blade 7 Lag Measurement 

678 

R ^7 

Noise covariance 

679 

b 57 

Bias 

680 

. k C7 

Scale factor error 

681 



682 

s 

Noise covariance 

683 

v 

Scale factor error 

684 

685 


Bias 








Table 3.4 (Continued) 




Generally, most of the parameters which are identified in 
a run are those used in forming the polynomial expansions of 
the aerodynamic coefficients. Parameters 1 through 400 are 
reserved for this purpose. However, the user should always 
use the low end of this range, when all locations are not 
needed, as this will reduce the run time. Chapter IV will 
explain in more detail how these expansions are accomplished. 




CHAPTER IV 
PROGRAM INPUT 


Two forms of input are used by the NLSCIDNT program. Infor- 
mation about the parameters and the selection of various options 
are read from cards. The time histories of the measurements and 
controls are read from a mass storage device, such as magnetic 
tape or disk. 

4.1 CARD INPUT 


All data cards are read by subroutine INPUT. The forms of the 
input fall into several types, which are listed in Table 4.1. 

The input sequence is as follows: , 

(1) One card of type 1. The information on this card is 
printed at the top of every page of printout. 

(2) One card of type 2. For an explanation of "outer" and 
"inner" iterations, see Section 5.2.2. A "prediction 
run" is one in which the trajectory is simulated using 
the control input time histories and input parameter 
values; no identification is attempted. 

(3) One card of type. 3. This card sets various option 
flags. If IPL0T = 2 or 3, the user must provide a 
tape on logical unit 3 (see Section 5.3). 

(4) Two or more cards of type 4. These cards tell the pro- 
gram which states are to be integrated from the equa- 
tions of motion, which states are to be looked up from 
an input time history, which measurements are used and 
their order, and which expansion variables are used and 
their order. If no states are to be looked-up, then a 
card beginning with "READ" is not necessary. Any state 
variable whose value does not appear on a "STATES" or 
"READ" card will automatically be held constant by the 
program at its initial value, which is one of the input 
parameters. These cards may appear in any order except 
the "*" card must be last. 


49 




Table 4.1 
Input Card Types 


■ 

COLUMNS 

FORMAT 

VARIABLE 

NAME 

DESCRIPTION 

1 

1-60 

15A4 

TITLE 

Program identification information. 

2 

in 

i 

12 

K1MAX 

Maximum number of outer iterations 

>0, normally 

= 0, prediction run 

=-l, only Oth iteration will be 
computed (usually selected 
for diagnostic purposes only) 

. 

9-10 

12 

K2MAX 

Maximum number of inner iterations 
> 0, normally 
= 0, prediction run 


12-15 

14 

NN 

Number of data sample points (< 501) 


16-25 

E10.3 

DELTA 

Time interval between sample points 
of the time history data (in sec.) 


26-35 

E 1 0.3 

RELER1 

Relative error bound for x and 

A A 

9x/90 allowed in integration algo- 
rithm (defaults to 10" 5 ). 


36-45 

E10.3 

i 

ABSER1 

Absolute error bound (defaults to 
1. 7 x 10** 4 ) . 


45-65 

E10.3 

FMARQ 

Initial value of Marquardt para- 
meter (Defaults to 1.) 

■ 

56-65 

E10.3 

XMARQ 

Factor by which the Marquardt para- 
meter will be increased or decreased 
as search proceeds. (Defaults to 2.) 


66-70 

15 

LLL 

=1, the Marquardt parameter is com- 
puted as the current value of 
FMARQ times the largest element 
of the gradient (9J/90). 

=0, the Marquardt parameter is 
FMARQ. 


SO 












Table 4.1 (Continued) 


BBSKi 

1 

COLUMNS 

FORMAT 

VARIABLE 

NAME 

DESCRIPTION 

3 

5 

11 

IPL0T 

=3, 

if both printer plots and a tape 
storing the time history data 
are to be made 





=2, 

if only the tape is to be made 





= 1, 

if only the printer plots are to 
be made 





=0, 

if neither 


10 

11 

IPLTC 

=1, 

if printer plots of control 
input time histories are to be 
made 





=0, 

if not 





(IPLTC is ignored if IPL0T = 0 or 2) 


12-15 

14 

I DATA 

=+k 

then the first k sample points 
of the measurement and control 
time histories read into the 
program will be printed 





=-k, 

then every kth sample point 
wi 1 1 be printed 





= 0 

if they are not to be printed 


20 

11 

IPRNT 

Sets level of detail of the diag- 
nostic printout (see Section 4.2) 





=3, 
=2, 
= 1, 
=0, 

highest level 
medium level 
lowest level 
normal printout only 


22-25 

14 

INCPR2 

=k. 

if x, 3x/30, y, 3y/30, and y are 
to be printed every kth sample 
point 


26-30 


INCPR3 

Default value is 50. INCPR3 is ig- 
nored if IPRNT < 3. 


33-35 

13 

INCPLT 

=k, 

every kth data point is plotted 
in the printer plots (default 
=1) 
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Table 4.1 (Continued) 


COLUMNS 


FORMAT 


VARIABLE 

NAME 


IRCMP 


IRVARY 



DESCRIPTION 


=1, a priori information matrix is 
to be read 

=0, a priori information matrix is 
not to be read but assumed zero 

=0, estimated noise covariance, R, 
to be computed as the diagonal 
of the innovations covariance 
(see Appendix A) 

=1, R computed from parameters 

=1, R is assumed nonstationary (it 
is a function of time which user 
must code into subroutine MEAS). 

=1, R is assumed stationary (con- 
stant throughout an iteration). 

Maximum number of integration steps 
permitted in one sample interval. 

If non-zero, it is the number of 
parameters to be identified on this 
run. 

Flag to use linearized state and 
measurement equations if f 0. 


= STATES, if the card lists the input 
names of states which are to be 
integrated 

= READ, if the card lists the input 
names of the states which will 
be looked-up from input time 
histories 

= MEASURE, if the card lists the 

input names of the measurement 
variables 

= EXP VARS, if the card lists the 
input names of the expansion 
variables 

= *b, if no more cards of this type 
are to be read ("b" indicates 
a blank). 

Note that only the first two charac- 
ters are actually read; these words 

must be left justified. 









Table 4.1 (Continued) 



— 

■ 

COLUMNS 

FORMAT 

VARIABLE 

NAME 

DESCRIPTION 

4 

11-80 

14(A3,2X) 

(LL(J), 

J-1.14) 

Input names of state, measurement, 
or expansion variables as listed in 
Table 4.2. Note that only the first 
three characters are actually read; 
these words must be left justified. 
Blank words may appear between input 
names . 

5 

1 

' 

A1 

ECHK 

= blank, if this card contains para- 
meter information 

= *, if no more cards of this type 
are to be read 


2-4 

13 

J1 

Parameter index (see Table 3.4) 


5 

A1 . 

J2 

= *, if this parameter is to be 
identified 

= blank, if not 


6-11 

A6 

PLABJ1 

Label for the parameter for printout 





purposes 


12-29 

i 

D18.0 

PJ1 

Initial value of the parameter 


31-32 

13 

. 

11(1) 

=k, this parameter appears in the 
polynomial expansion of the kth 
aero coefficient (see Table 3.2), 

Ignored if J1 > 400. 


34-38 

511 

11(0), 
J=2 , . . . ,6 

Exponents of the expansion variables 
in the aero coefficient polynomial 
term containing this parameter. 


51-60 

D10.0 

PLJ1 

Lower bound of the parameter. 


61-70 

D10.0 

PUJ1 

Upper bound of the parameter. (If 
both PLJ1 = 0 and PUJ1 =0, then. 

PLJ1 defaults to -10 27 and PUJ1 de- 
faults to +10 27 ). 
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Table 4.1 (Continued) 


COLUMNS 

FORMAT 

VARIABLE 

NAME 

DESCRIPTION 

1-80 

8F10.0 

INF0(J,K) 

K~ X y • • • jID 

Jth row of the a priori information 
matrix; continue on another card if 
necessary (m = number of parameters 
being identified). 

1 

A1 

JA 

= T, the label which follows on this 
card is for time (the indepen- 
dent variable) 

= Y, the label is for a measurement 
variable 

= U, the label is for a control 
variable 

= X, the label is for a state vari- 
able (only look-up states are 
affected) 

3-4 

12 

. 0 

The index number of the variable. 
If JA = T, then J is ignored. 

7-30 

6A4 

NAME 

A label of six 4-character words. 
Usually the variable name appears in 
the first 3 words and its units in 
the 1 ast 3 words . 












(5) As many cards of type 5 as needed to input all informa- 
tion about the parameters; each card describes one para- 
meter. If a parameter is zero and is not identified, 

no card need be inserted for it. A card with in 

column 1 tells the program that all cards of this type 
have been read. The system of units for the parameter 
values is feet, pounds, slugs, seconds, and radians 
or meters, newtons, kilograms, seconds, and radians. 

It is the responsibility of the user to maintain con- 
sistent units throughout. 

(6) Cards of type 6 are read only if IINF0 = 1 on card type 

3. These cards are used to read the a priori information 
matrix. This matrix must be mxm, where m is the 
number of parameters being identified. Furthermore, 
the order of the rows and columns must correspond to the 
order of the identified parameters as listed in the 
cards of type 5. The matrix is read by rows; each row 
begins on a new card and continues on more cards as 
needed. 

(7) As many cards of- type 7 as needed to alter the labels 
of time, measurements, and controls in the printout and 
printer plots of the time histories. These labels 
default to the output labels listed in Table 4.2. There- 
fore, the most common use of these cards is to change the 
units shown in the labels. If none of the labels need 
altering, then no . cards of type 7 are necessary. 

4. 2 EXAMPLE INPUT DECK . 


As an illustration of the card input, consider the following 
example. The data processed were simulated from a rotorcraft in 
level flight at 100 knots at sea level . Longitudinal cyclic pitch 
control inputs excited the longitudinal dynamics. However, some 
perturbations in the lateral dynamics were observed and of enough 
significance that they could not be neglected. A NLSCIDNT run 
was made to identify three aerodynamic derivatives. 

Figure 4.1 shows the input card deck for this example run. 

The first card contains the title of the run. The second card 
specifies one outer iteration, a maximum of four inner iterations 
permitted, that there are 101 sample points, and that 0.02 second 
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Table 4.2 

Input Names and Output Labels of Program Variables 



SYMBOL 

mimwrn 

1 

OUTPUT 

LABEL 2 

TIME 

t 

— 

TIME 

(SECONDS) 


u 

u 

LONGIT VEL 

(M/SEC) 


V 

V 

LATERAL VEL 

(M/SEC) 


w 

w 

VERTICAL VEL (M/SEC) 


p 

p 

ROLL RATE 

(RAD/SEC) 


q 


PITCH RATE 

(RAD/SEC) 


r 

R 

YAW RATE 

(RAD/SEC) 

c 

4> 

PHI 

ROLL ANGLE 

(RAD) 

J 

T 

9 

THET 

PITCH ANGLE 

(RAD) 

I 

A 

$ 

PSI 

YAW ANGLE 

(RAD) 

A 

T 

• 

B o 

• 

BOD 

BETA D0T 

(RAD/SEC) 

E 

B 1c 

• 

BCD 

BETA 1c D0T (RAD/SEC) 

c 

B 1s 

BSD 

BETA Is D0T (RAD/SEC) 


B o • 

BO 

BETA 0 

(RAD) 


B 1c 

BC 

BETA 1c 

(RAD) 


?ls 

BS . 

BETA Is 

(RAD) 



ZOD 

ZETA 0 D0T 

(RAD) 


*lc 

ZCD 

ZETA 1c D0T (RAD/SEC) 


*1. 

ZSD 

ZETA Is DOT (RAD/SEC) 


“s 

OME 

OMEGA PERT. 

(RAD/SEC) 



10 

ZETA 0 

(RAD) 


? lc 

ZC 

ZETA 1c 

(RAD) 


^ls 

ZS 

ZETA Is 

(RAD) 



RSP 

PSI-S PERT. 

(RAD) 


56 


















Table 4.2 (Continued) 




























Table 4.2 (Continued) 


SYMBOL 



OUTPUT LABEL 2 


FLAP ANGLE-BLADE 4 (RAD) 
FLAP ANGLE-BLADE 5 (RAD) 
FLAP ANGLE-BLADE 6 (RAD) 
FLAP ANGLE-BLADE 7 (RAD) 


LAG ANGLE-BLADE 1 (RAD) 
LAG ANGLE-BLADE 2 (RAD) 
LAG ANGLE-BLADE 3 (RAD) 
LAG ANGLE-BLADE 4 (RAD) 
LAG ANGLE-BLADE 5 (RAD) 
LAG ANGLE-BLADE 6 (RAD) 
LAG ANGLE-BLADE 7 (RAD) 


(cos i|> R ) m COS 

(sin Vm SIN 


COS ROTOR AZIMUTH (RAD) 
SIN ROTOR AZIMUTH (RAD) 









Table 4.2 (Continued) 



SYMBOL 

B!H m 

mmm 

OUTPUT LABEL 2 

c 

9 o 

TO 

COLL PITCH (RADIANS) 

0 

8 lc 

TC 

LAT CY PITCH (RADIANS) 

N 

8 ls 

TS 

LON CY PITCH (RAD IANS) 

T 

R 

6 tr 

DTR 

TAIL COLL (RADIANS) 

0 

5 

e 

DE 

ELEVATOR (RADIANS) 

L 

6 

a 

DA 

AILERON (RADIANS) 

S 

' 6 

r 

DR 

RUDDER (RADIANS) 


6 f 

DF 

FLAPERON (RADIANS) 


a 


ALPHA 

E V 

B 

BETA 

BETA 

X A 

P 

P 

ROLL RATE 

P R 
A I 

q 

Q_ 

PITCH RATE 

N A 

r 

R 

YAW RATE 

SB 

e o 

TO 

COLLECTIVE 

I L 
0 E 

e lc 

TC 

LAT CYCLIC 

N S 

e ls 


LONG CYCLIC 


5 tr 


TAIL ROTOR 


5 e 

DE 

ELEVATOR 


V 

DA 

AILERON 


S r 

DR 

RUDDER . 


S f 

DF 

FLAPERON 
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Table 4.2 (Concluded) 



SYMBOL 

1 1 
iqih 

OUTPUT LABEL 2 

■■■ 

V 

MU 

MU 


u 

u 

u (NORM) 

V 

V 

v (NORM) 


w 

w 

w (NORM) 

■ 






*The underlined characters of the input name are 
what the program actually reads. 

2 "MS" in the measurement labels indicates the 
measured quantity; "ES" appears instead in the 
label to indicate the estimated quantity. 
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COMMENTS 


CARD CARD 
TYPE NO. 


DATA CARDS 


1 

2 

3 


6 

7 


i. 

1 . 

o. 


| ATffUl. rul'lSF f>VM*«ICS fc*AMP|f.«- 

I it t o i n~o ? -4 

• i II 3 25 2S 

STAIEjj, V P R PH! 


0 


5. 

Of AHJ 

w o Theta 


— 

6,~ 

mE A SURF t # 

ay beta p r 

PHI 

4 — — 


7. 

MP VApfii 

ftf T A OSP 

■ ■- 

- 

< 

♦ 


J 

f 

9. 

i r Y 

°. 

6 


10, 

? CYV 

-*onq?6ft 

7 


l 1 . 

3 rvp 

- .ClftftQ? 

ft 



« cvp 

o . a « a 7 

9 


1 3 . 

S^rvo* 

.AollPl 

io ’ 


i a. 

ft C v O r 

O.I«6A 

1 1 


is.. 

7 ri 

n . 

i7‘ 



p nv 

- .noosJSO., 

i« 


• 7 * 

9*C| P 

•• 0 . ss 

1 *» 

.10000000* o • 

1*. 

in n r 

0, t Sui 

?0 



11 f 1 OA 

0.1 *S> 

?1 


?o . 

I? ri or 

• 0?7|)t 

?? 



i 3 r wo 

0. 

?ft 



1 <i ruv 

•00 oss 3 

?9^ 


?3. 

IS r mp 

-.061S3 

30 



to rt.p 

*0.1 H 03 

31 


?S. 

17 tun a 

•0OS9P0 

3? 



10 tMI'O 

-o.n?i 

33 


P 7 . 

1 9 fh03 


?6 30^ 



2o#ri n$ 

0. 

|7 oi 



?i *n nsl 

0. 

1 7 0| 



ft if PHO 

•00P3769 




ft? TRAP 

6.S 



3?. 

B3 ft 

6?0. 



33* 

ftft R 

6S. 



■3 fl. 

ftS *'ASft 

I'l? . 



3S. 

BO I* 

1 6 « <> (1 A 



3ft. 

ft 7 fy 

?'t9no t 



37, 

ftfl 17 

J7SS0. 



3*. 

ft 9 1*7 

0. 



39. 


1117. 



00. 

101 l»0 

179.7* 




103 Vi o 

*».«2 




t 

] 

1 


«3. 

7 70?. « 

?nl*>.l 5?.f.5i | 

[ 



?0|9. i 

692l.il. 176,60 j 

[ * 


OS. 


171.. 6 0.07228 J 

1 . 


00. 

T T t M f 

fsrci rff in iso 

•A 4- 


07. 

7) pHf 

( P AO ) 

r 



. ti tie 

max iterations, max step cuts, no. of 
sample points, and sample interval 
option flags 

specification of variables used in 
this run 


-parameter number 
-parameter name 
-parameter value 
■upper bound 
-lower bound 

index of coefficient in whose expansion 
this parameter appears 

exponents of expansion variables 


indicates the parameter is to be 
identified 


a priori information matrix 
specification of output labels 


Figure 4.1 Example Input Card Deck 


is the sample interval. The third card asks for printer plots, 
for the first 11 sample points of the input time histories to 
be printed, for the highest (third) level of diagnostic print- 
out at intervals of 25 sample points, and for the second level 
of diagnostic printout at intervals of 25 sample points. It 
also specifies that an a priori information matrix is not to be 
read, that the noise covariance R be computed fiom the para- 
meters, and that R is constant with respect to time. The 
maximum allowable number of integration steps per sample 
interval is set to default to SO. 

The fourth card tells the program which states are to be 
integrated. The states which do not appear in the fourth card 
will be held constant at their initial values. The fifth card 
tells which states are to be read in from the test data. 

The sixth card specifies the measurements available to the 
run and their order. The seventh card indicates the expansion 
variables. 

The first 21 parameters (specified on cards 9 through 42) 
are used in the expansions of the aero coefficients which appear 
in the equations of motion. 

The parameters which are to be identified are P(9), P(20), 
and P(21). P(16) is bounded such that the identified value 

will not be permitted to become positive. 

Cards 43 to 45 contain the a priori information matrix. 
Cards 46 to 47 give new labels to time and 4> . 



MASS STORAGE INPUT 


4. 3 

The time histories of the measurements, controls, and states 
which are looked-up rather than integrated or fixed are read by 
subroutine INREAD. The subroutine INREAD supplied with the pro- 
gram assumes these data are stored on a mass storage device (tape, 
disk, etc.), referenced by F0RTRAN logical unit 2 and read by an 
unformatted READ statement: 

D0 10 K= 1 , NN 

10 READ ( 2 , END=900) T (K) , (Y ( J , K) , J®1 ,NP) , (U ( J , K) , J=1 ,NQR) 

where 

NN = number of sample points, 

T(K) = the time at the K-th sample point, 

Y(J,K) = value of the J-th measurement at the K-th 
sample point , 

• NP = number of measurement variables used in 

the run, 

U (J , K) = value of the J-th control variable or a 
state variable at the K-th sample time, 

NQ'r = number of control variables (4) plus the 

number of look-up states (those not inte- 
grated or fixed) 

T, U, and Y are single precision arrays. 


The user must ensure that the time, history data conform to 
these rules: 


(1) The order of the measurements must be the same as their 
specified order in the card input (card type 4 in 
Table 4.1). 

(2) The controls must always be in the order of: 


• collective pitch 

• lateral cyclic pitch 

• longitudinal cyclic pitch 

• tail rotor pitch 


• elevator 

• aileron 

• rudder 

• spoiler 


(as they appear in Table 3.1); each must be included, 
even if it is always zero. 


(3) Look-up state variables must be appended as if they 
were additional controls. Their order must be the 
same as their specified order in the card input (card 
type 4 in Table 4.1). 


(4) No more than four look-up states are permitted. 

For the example run described in Section 4.2, the k-th logi- 
cal record of the time history data must be: 

t k * a Y^ t k' ) ’ 3 m (t k )> ’ 5 e ^k' ’ 


6 a^H' }> V^k^’ 5 sp^kV» q( ^ t k ) ’ 8 (t k^ 


The user will almost certainly have to rearrange his data to 
conform to the structure required by this subroutine. If he 
believes he will process the same data repeatedly, it would be 
worthwhile to perform this rearrangement once before the first 
NLSCIDNT run and take advantage of the efficiency of this INREAD. 
However, if the data, will be processed only a very few times, he 
may choose -instead to write his own subroutine INREAD, which 
will accept the data in their original form. Figure 4.2 shows 
the source code of the INREAD subroutine supplied with the pro- 
gram, which may serve as a model for the user's own. 


sumjouMhf inhfao (u»no,y,nhi r.uN.NOio 

niK-LNSlUN U(nUiHN) »Y<M‘.UN) »T(NN) 

COXMOh /CUl-IOll/ 10 1 L I , J NT A*'t 

00 10 K=l,HH 

Of Al>< |MA|*f .E.HR=90n,FNl>s90?) T (K ) » ( Y C J . K ) , J= 1 1 HP) » (UCJ i K ) . J = I . NOB) 
10 C(ir:l iitUt 

99 IHtUHiJ 

900 kill (I- (| O.i 90 U K.T(K-l) 

901 » ull.tA |C l»M. • ♦♦♦ ERROR 1 M SUHHOUT I Nt InRE'AD, DEVICE ERROR CtltCTEO 
1 *hll.f aTIEmIMINC 10 HF A 1 ' OATa POINT I .Hi • . '/fiXf I IHE TIRE. AT THE 
?PH»vlHUS* l*'JHI» I(A-I), k aS ' i I R E 1 0 , 6 ) 

5 I OH 

90? •**! |F (1.0.903) KfT(K-l) 

901 FilR^Al ( I >10. i *« * ERROR IN SU|,R(HJTlNf 1 nH£aD. tND OF FILE ENCOUNTER 
ICO r.H|l.F a 1 1 1 kR 1 1 NC TCI read DaIA POINT • . H • ' . ' /5X . I THE TIRE AT Th 
2E PREVIOUS RO|Nf, l(K-ni WAS IiIPEH.6) 

5 I 01* 

FnO 


Figure 4.2 Subroutine INlUiAD Source Code 




V. PROGRAM OUTPUT 


Most of the program's output is written by the printer, in- 
cluding optional plots. The program also has an option for writing 
data on a mass storage device, which gives the user the information 
necessary to make off-line, pen-and-ink plots (e.g. Calcomp plots). 

5.1 PRINTED OUTPUT 

The information printed by the program falls into six classi- 
fications. Not all classifications may appear since the user has 
considerable control over the detail of the printout through his 
selection of various options on input card type 3. The six classes 
are : 

(1) repetition of the inputs to the program, 

(2) information on the integration of the differential 
equations , 

(3) information on the Levenberg-Marquardt iterative search, 

(4) information on the steps taken in the identified para- 
meters and their error standard deviations, 

(5) the final parameter estimates (always given) , and 

(6) printer plots showing: (a) the fit of the estimated 

measurements to the actual measurements (IPLOT = 1 or 3), 
and '(b) the control time histories (IPLTC = 1) . 

5.1.1 Example of Printed Output 

Figure 5.1 consists of selected pages from the printout of 
the same example run whose input was discussed in Section 4.2. 

It contains examples of all the available printed output. 

The information in Figure 5.1a is a repetition of the input 
and is self-explanatory except for the flag arrays. .ISFLAG(J) 
contains: 0 if the jth state variable will be fixed to its initial 
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Figure 5.1a Example NLSCIDNT 


Printout 


value, +k to indicate it is the kth element in the vector of 
state variables to be integrated, and -k to indicate it is the 
kth element in the vector of look-up states. NS FLAG contains 
the indices of the states to be integrated. NMFLAG contains the 
indices of the measurement variables used in the run. NZFLAG 
contains the indices of the expansion variables used in the run. 
NZFLAG (1) will default to 1 if no expansion variables are used. 

Figure 5.1(b) shows the parameter input information. If a 
parameter does not appear in the list, its value is zero and it is 
not identified. Figure 5.1(c) presents a sample of the input 
time history data. Only the first 11 sample points were printed 
because IDATA was set to 11 by the input deck (refer to card type 3 
in Table 4 . 1 and Figure 4 . 1) . 

The CV matrix in Figure 5. Id is. 

' N 

CV = s 
i-i 

where 

y^ = measurement vector at t * t- , 

u i = control vector at t = t^, 

x Li = vector of look-up states at t = t^, 

N = total number of sample points. 

Eleven elements are zero because all controls except longitudinal 
cyclic pitch are zero for all N points and because the four 
look-up states are not used. 

The IU array contains the indices of non- zero control 
variables. Look-up states would be treated as additional controls. 
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Figure 5.1b 
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The initial noise covariance matrix, R, is also printed in 
Figure 5.1(d). Because none of the instrument noise covariances 
were specified by the input deck, each of the diagonal elements 
of. R was approximated as 



0.05 

N 


N 

X 

i=l 


C/jCti)] 


(5.2) 


That is, the variance of the noise in the jth instrument is 
approximated as 5% of the variance in the instrument's signal. 


. Information on the integration of the state and sensitivity 
differential equations for the zero-th iteration of the Newton- 
Raphson search begins in Figure 5.1e. A "matrix" of initial con- 
ditions appears in this figure. The first column is the initial 
condition of the x vector (estimated states) . This column has 
length four because only four of the nine states. are being inte- 
grated. , The second column is the initial condition of ax/30^, 
i.e. , the sensitivity of the estimated states to the first identi- 
fied parameter, which is P(9) in the example. The third column 
is the initial condition of 3x/S6 2 > etc. The second and succeed- 
ing columns should always be zero unless the initial conditions 
of the states are among the identified parameters. 


The four groups of values printed between the first line of 
dashes and the line of stars in Figure 5.1e are intermediate 
results in computing the time derivatives of x and 3x/30. First 
time (T) and the expansion variables (2) at that time are printed. 
In the example, there are two expansion variables. Next is the 
u-vector showing the controls and look-up states (appended as addi- 
tional controls) at time T. The four normal controls are zero. 

The first look-up state, x(5) = w, equals 9.4202. The second 
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Figure 5.1e 


&nd third look*up stages, x(5) = q and x(8) =9, are zero. A 
potential fourth look-up state is printed, even though it is not 
defined. 

The x-vector printed next is the nine-element state vector. 
For example, x(l) s u - 179.75 and x(3) = w = 9.4202. Next are 
printed a series of terms which appear in the equations of motion. 

V W 

BVP = P ' 

co 

CBVQ * -J- q 
• ■ * 00 

■QD-isVi-q. ■ 


QB = q^sb 

QC - q^Sc 

BVR = 7V^ r (5 • 3 

The 33 aerodynamic coefficients .(C) are printed in order across 
the page; hence, C(ll) = Cy^ r =0.3468. Last in this group appear 
the time derivatives of the state variables being integrated, DXI . 
In this example, there are four: v, p, f, and 6. 

The second group contains information on sensitivities with 
respect to the first identified parameter, which is P(9) = . 

Time is repeated and followed by 2TH , the partial derivatives of 
expansion variables with respect to P ( 9 ) . Next - are the partial 
derivatives of the quantities defined in Eqs. (5.3). The partial 
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derivatives of the aero coefficients are "DC." Note that 
DC (19) a 1 because C(19) » Cj^ = P C 9 ) - Last in this group is 
DXTHI , which is d [ 3x/3P (9) ] /dt for the four integrated state 
variables. 

Groups three and four are similar to group two except that 
all partial derivatives are with respect to P(20) and P(21), 
the second and third identified parameters , respectively . Note 
specifically that at T * 0, DC (17) » 0 in both groups three and 
four, even though 

c Cl") ■ P(7) * P(20)(i,-2 2o ) 1 * P(21) (2 2 -: 2o ) 3 

- PC?) * PC20)o sp ♦ P(21)«j p (5.45 

This is because the partial derivative of C(17) with respect to 
either P(20) or P(21) is zero when = 0, which it does 

here. 

The information printed between the line of stars and the 
second line of dashes in Figure 5.1e shows the results of the in- 
tegration at the j-th sample point. In this instance, j = 1. 

In the middle of the line of stars in Figure 5.1e is printed 
KOUNT, which is a running total of the number of times subroutine 
STEP has been called. This is approximately the number of times 
the states' and sensitivities' time derivatives have been computed 
More will be said of its importance in Chapter VI. 

Next is ”XHAT(1)" which is the vector of integrated states 
and their sensitivities with respect to the identified parameters 
at the (in this case) first sample point (t^ * 0). They appear 
in the order 


x l» x 2 * *3’ x 4 ’ 3x 2 / 3 6 1 , ax-/36 1 , 3x 4 /36 1 , 


3x^/ c9 7 , ..., Sx^/36. 


76 



They are followed by "YHAT” , the measurement estimates and their 
sensitivities with respect to the identified parameters in similar 
order. M Y- VECTOR" is the vector of actual measurements at the 
same time. 

As the integration proceeds, the groups shown in Figure 5.le 
are repeated at intervals specified by the program's input. These 
pages are omitted from the figure to save space. Figure 5. If has 
the very end of this stream of information at its top. 

The RCMP matrix is printed next.. Its diagonal elements are 
computed as 


RCMP . . 
1,1 


_ 1 


N 


n * ^V-W 1 


(5.5) 


These are the variances of the errors between the actual and the 
estimated measurements; and, therefore, are a measure of the good- 
ness of the fit. Also, it is the best estimate- of the measurement 
noise covariance (see Appendix A] . 

. Figure 5. lg presents the results of the zero-th (i.e., initial) 
iteration of the Newton- Raphson search. The DJ vector is the 
partial of J, the negative log- likelihood cost function, with 
respect to the three identified parameters. The D2J matrix is 
the 3 J/39j36^ matrix. Also printed are its inverse,- its eigen- 
values, and its eigenvectors. The significance of the number of 
"eigenvalues effecting a full-sized step" will be deferred to 
Chapter VI. 

Next are printed the values of the parameter estimates, the 
lower bounds on the error standard deviations of these estimated 
values, and their F values, which are the ratios of the squares 
of the parameter values to the squares of their respective standard 
deviations. The value of the likelihood function and the length 
Ox the gradient are printed. Of course, it is the likelihood func- 
tion that the program is seeking to minimize, and the length of 
the gradient will be zero at a minimum. 
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^ top of Figure 5 . 1 h are the "new parameter values" 
found by adding the "step" to the "parameter values" printed 
ea,r ^ 1 ® r * These new parameter Values are presumed to result in a 
value of the likelihood function less than the one obtained for 
this iteration, and they are the parameter values which will be 
used at the beginning of the next iteration. "Normalized step" 
is the ratio of the step vector to the parameter values in Figure 
5 . 1 g, except that division by a zero parameter value is replaced by 
division by one. This ratio tells at a glance the relative size 
of the step, and it will decrease as a minimum is approached. 

R is the measurement noise covariance which will be used 
in the next iteration. What follows it is actually part of the 
next iteration and parallels the printout of Figure 5.1e 

This example run was allowed only one iteration and, there- 
fore, the parameter search aid not converge. Thus, the message 
"maximum number of iterations exceeded" appears in Figure 5.ij. 

The matrix A is the sum of the a priori information matrix, 

M ao’ and the information matrix from the last iteration. M 
The vector B -is the sum "" c 


B = M P + M, P 

a P o last last 


( 5 . 6 ) 


where 


Finally, 


last 


the initial (input) values of the identified 
parameters, and 

the values of these parameters used in the last 
iteration (not the last "new parameter values") 


PSTAR — A B ® (M +M.. ) ^ (M p + m P 1 ft 7^ 

/• ap last J u ‘ap r o ‘‘last^lasf 1 ( 5 * 7 J 
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which is the vector of the optimal parameter estimates which com- 
bine the information obtained from this run and the information 
known a priori. Of course, to be valid this run should have 
been allowed to converge. A, B, and PSTAR will not be printed 
if the a priori information matrix is zero. 

Next are printed plots of the actual and estimated measure- 
ment time histories and the control time histories. Figure .5.1k 
is an example of the measurement plots and Figure 5.1$, of the 
control plots. In the measurement plots, "A" is always printed 
for the actual measurement and "B" for the estimated measurement. 
Control variables are plotted two to a graph; and, again, look-up 
states are treated as additional controls. 

5.1.2 Printout Control 


In most cases, the detail of the printout is controlled by 
the flag IPRNT in card type 3 of Table 4.1. For normal production 
runs this flag is set to zero, and the printout is shortened to 
essential results of each Newton-Raphson iteration. Setting 
IPRNT = 1, 2, or 3 results in increasingly greater detail of 
diagnostic information, as explained in Table 5.1. 

Setting IPRNT = 1 results mainly in added information regard- 
ing the computation of the parameter step in each iteration. 

Setting IPRNT = 2 results mainly in the printing of x, 3x/36, 

A 

y, and 3y/30 as the trajectory is integrated for each iteration. 
Finally, setting IPRNT = 3 results in detailed information on the 
trajectory integration. 

5.2 MASS STORAGE OUTPUT 

An option is provided to write on mass storage: (a) the last 

vector of identified parameters, (b) the last information matrix, 
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Table 5.1 

Printed Information Controlled by IPRNT Flag 


FLAG VALUE 

INFORMATION PRINTED 

COMMENTS 

IPRNT = 0 

1. Card input data repeated 

Only parameters mentioned 
in the input deck are 
printed under "Initial 
Parameter Values." The 
a priori information matrix 
is printed only if it is 
non -zero. 


2. For each iteration of the Newton- 
Raphson search: 



a. Parameter values at start of 
iteration 

b. Parameters' standard devia- 
tions 

c. Parameters' F-values 

d. Value of negative log likeli- 
hood function, J 

e. Length of first gradient of J 

f. - New parameter estimates 

g. Parameter steps 

h. Innovation covariance (RCMP) 



3. Final combined parameter set 
(PSTAR) 

| Only if a priori informa- 
tion is non-zero 

IPRNT = 1 

Same as for IPRNT =0 plus: 



1. Values in ISFLAG, NSFLAG, NMFLAG, 
NZFLAG, and IU arrays 

These arrays described in 
Section 5.1.1 


2. CV matrix 

3. Initial R matrix 

Described in Section 5.1.1 


4. For each iteration of the Newton- 
Rapnson search: 



a. First gradient of J (DO) 

b. Second gradient of J (D2J) 

c. Eigenvalues and eigenvectors 

of (D2J)”* and number of 
eigenvalues effecting a full- 
sized step 

Significance described 
in Chapter VI 


5. Combined information matrix (A) 

Only if a priori informa- 
tion matrix is non- zero. 
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Table 5.1 (Continued) 


FLAG VALUE 

INFORMATION PRINTED 

COMMENTS 

IPRNT a 2 

Same as for IPRNT * 1 plus for each 
iteration of the Newton-Raphson 
search: 

a. States and their sensitivity esti- 
mates (XHAT) 

b. Measurements and their sensitivity 
estimates (YHAT) 

c. Actual measurements (Y-VECTOR) 

d. Number of calls to subroutines 
STEP, C0EF, AND DC0EF 

e. Inverse of D2J matrix 

) Printed for the first 
( 4 sample points and for 
( every (INCPR2)-th 
J sample point 

Significance described in 
Section VI. 

IPRNT = 3 

Same as for IPRNT -2 plus for each 

iteration of the Newton-Raphson 

search: 

a. Initial conditions of states and 
their sensitivities 

b. Time and values of expansion 
variables (T and Z) 

c. Control vector including look-up 
states (U- VECTOR) 

d. Nine-element state vector 
(X- VECTOR) 

e. Values of V, BVP, CBVQ, QD, QS, 
QC, BVP, BVR 

f. Aero coefficients (C) 

g. Time derivatives of integrated 
states (DXI) 

h. Partial derivatives of variables 
in b through g above with respect 
to the identified parameters 

I Printed for the first 
V 3 sample points and for 
/ every (INCPR3)-th 
1 sample point 
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(c) the time histories of the actual and estimated measurements, 
and (d) the time histories of the control variables. The mass 
storage device (tape drive or disk) is currently designated as 
logical unit 3. The data are written in unformatted records by 
statements equivalent to the following: 

WRITE (3) NS , NQ , NP , NN , M, (PID(J) ,J=1,M) , ((INF0LCJ.K) , 

1 J=1,M) ,K=1,M) 

D0 10 K= 1 , NN 

10 WRITE (3) T(K),(UCJ,K),J=1,NQ),CY(J,K),J=1,NP), 

1 (THAT (J,K) ,J=1,NP) 

ENDFILE 3 


where 


NS 

NQ 

NP 

NN 

M 

PID 

INF0L 

T 

U 

Y 


number of states integrated in the run* 

number of control variables in the model plus 
number of possible look-up states (currentlv 
NQ=8) ' 

number of measurements used in the run, 

number of sample points in the run, 

number of parameters identified in the run, 

vector of identified parameters used in the 
last iteration of the run, 

last information matrix, 

array of time values at sample points, 

array of control and look-up state variables 
at sample points , 

array of actual measurements at sample Doints 
and ’ 


YHAT * array of estimated measurements at sample points. 
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CHAPTER VI 


EFFECTIVE USE OF THE NLSCIDNT PROGRAM 


This chapter contains information on selected topics and 
guidelines for effective use of the NLSCIDNT Program. 

6.1 INTEGRATION ALGORITHM 

Each outer iteration of the Levenberg-Marquardt search requires the 
integration of n(m+l) ordinary differential equations, where 
n is the number of integrated states and m is the number of 
parameters identified. This is the most time-consuming task the 
program must perform. Most of this time is spent evaluating 
dx/dt and d(3x/30)/dt. Therefore, an integration algorithm 
which requires relatively few evaluations of these derivatives is 
more efficient than one requiring more. 

Adams formulas are well known for this desirable character- 
istic, but most implementations are not self-starting (for ex- 
ample, four Runge-Kutta steps may have to be taken first), they 
have a fixed order (usually they require four previous points), 
and they have clumsy machinery for varying the stepsize (which 
may cost them much of their inherent advantage in speed). A 
recent implementation of an Adams formulation due to Shampine and 
Gordon [3] overcomes all these objections. It is called a PECE 
method because for each step it: 

• predicts- the solution of the differential equations at 
the end of the step using the Adams -Bashforth predictor 
of order k, 

• evaluates the derivative at this predicted solution 
point, 

• corrects its solution value at the end of the step using 
the Adams -Moulton corrector of order k+1, and 

• evaluates the derivative again using the corrected 
solution. 
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The PECE method was designed to be self -starting. Also, 
the order and the stepsize are both varied routinely to reduce 
the number of evaluations of the derivatives while meeting relative 
and absolute error bounds imposed on the solution by the user. 

These bounds are imposed as follows. Let y^ be any element 
of the x and 3x/36 vectors as computed by the PECE method 
at the' i-th point, and let be the difference between y^ 

and its true value, y^, at that point. Because y^ is unknown, 
can only be estimated. However, the theory says that n^, 
the difference between the predicted and the corrected values for 

a 

y^ , is an upper bound on the error e^. That is, 

I | < | n i | (6.1) 

The program bounds the error by requiring 

I n i | < RELERRMyJ + ABSERR (6.2) 

where RELERR and ABSERR are inputs to the program. If the user 

- 5 - 4 

does not supply values, they default to 10 and 1.7x10 , 

- 4 

respectively. (The value 1.7x10 is 0.01 degrees expressed in 
radians.) Note that choosing ABSERR=0 will cause difficulties if 

A 

the solution y^ passes through zero. 

If the error bounds are unnecessarily restrictive, then 
the program will be forced to take many small steps and waste 
time. If the error bounds are too lax, the solution of the dif- 
ferential equations will be inaccurate, and the subsequently 
computed step in the identified parameters will also be in error. 
The user should experiment to find values of RELERR and ABSERR 
which are best suited to the problems he solves. 

The program constrains the order of the predictor equation 
between 1 and 12. That is, as few as one and as many as twelve 
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past time points of data may be used to predict the next. Obvi- 
ously, this requires setting aside sufficient storage for the 
maximum of twelve past points, which is 12n(m+l) 
locations. Thus, higher orders quickly become impractical and pro 
duce diminishing returns anyway. In practice, twelth-order is 
seldom reached. 

The maximum permitted integration step-size is ten times 
the data sample step-size. The minimum is 0.02 times the sample 
step-size, unless another factor is read into the program (MAXNM1 
in Table 4.1). In general, the ratio of the integration step-size 
to the sample step-size or the inverse ratio is not an integer. 
Hence, the "mesh” points of the integration (i.e., the time 
points which mark the end of the integration steps taken) almost 
never coincide with data sample times. The values of x and 
Ox/30 must be found by interpolation between their values at 
mesh points. This is an inexpensive task, however, since the 
interpolating polynomial is readily available from the correction 
equation. Furthermore, the theory says that this interpolating 
polynomial is as accurate between mesh points as at them. 

In a typical run, the integration step-size is initially 
much smaller than the sample step-size. In fact, it may be 
repeatedly halved and the time derivatives of x and 3x/3e 
repeatedly computed until a small enough first step is taken to 
satisfy the error bounds. Once started, however, the integration 
step-size is steadily increased until it is 3 to 7 times the 
sample step size. Consequently, the run time is not directly 
proportional to the number of data points (see Figure 6.1). 

In the example run whose printed output is shown in Figure 
5.1, the initial number of integration steps per sample step was 
6.0. By the 100th sample point, this ratio was down to 0.52. 

The total number of steps taken is the number of times subroutine 
STEP is called. Because STEP may occasionally make more than one 
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Figure 6.1 Variation of Integration Step-Size in a 
Typical Run 


attempt at a successful step in order to satisfy the error bounds, 
the actual total number of times dx/dt and d(3x/38)/dt are 
evaluated may be slightly higher than twice the number of steps. 
(If STEP made no "false starts," it would be exactly twice.) 

The number of times subroutine C0EF is called reflects the number 
of evaluations of dx/dt. Because subroutine DC0EF is called 
once. in evaluating d(3x/30^)/dt for each 0^, i=l,...,m, it 

will be called m times more often than C0EF. 

The reader may now appreciate the efficiency of the PECE 
algorithm. In 101 sample points, it needed to evaluate the time 
derivatives only 118 times. The second-order Runge-Kutta algo- 
rithm would have required 202 evaluations of the time derivatives, 
for example. Greater savings are realized for longer data 
records . 

The cause and possible solution to error conditions in the use 
of subroutine STEP will be discussed briefly in Chapter VII For 
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more detailed information on the source of error conditions as well 
as more details on the theory and use of the PECE algorithm, the 
reader is referred to Ref. 3. 

6.2 CONVERGING TO A MINIMUM (AND OTHERWISE) 

6.2.1 Convergence Criteria 


After each Levenberg-Marquardt iteration, the program performs 
two tests to determine if it has found a minimum of the cost func- 
tion. If either test gives a positive result, the search ter- 
minates. These tests are: 

< 10 " 5 ( 6 . 3 ) 

is the cost after the current iteration and 
e cost after the previous outer iteration. 


(2) m V (ae) N (48) N < 10 ’ 3 (6-4) 

where (A3)^ is the "normalized step" vector (explained 
in Section 5.1.1), and m is the number of identified 
parameters. 

6.2.2 Step Cutting Procedure 

These criteria have served the purpose well in all cases 
except one- -step cuts. If the value of the cost function for an 
iteration is greater than the cost at the previous iteration, the 
program cuts back on the parameter step it has taken until a lower 
cost is found. Usually, this works satisfactorily, and the program 


( 1 ) 


J k ' J k-1 


where J, 


J k-1 is 


t 
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proceeds. However, sometimes round-off or other inaccuracies 
(such as error bounds which are too lax) cause the step to have 
been computed erroneously in the first place, and no amount of 
step cutting will result in a lower cost. The result is that the 
step is cut so much that one of the convergence criteria is met 
even though the cost may be nowhere near a minimum. The user 
should suspect this condition whenever he sees a large number 
of step cuts occurring just before the search appears to converge. 
A suggested remedy is given in Chapter VII. 

Nonetheless, the step cutting procedure in NLSCIDNT has 
proved itself very effective. It is a common occurrence that the 
cost function may be nearly insensitive to one or more linear 
combinations of the identified parameters. In such cases, the 
parameter step has large components in these directions of the 
parameter space. In Appendix A it is shown that the parameter 
step vector is the sum of vectors whose directions are the same as 
the eigenvectors of (3 2 J/39 2 ) ' 1 , i.e. , the inverse of the infor- 
mation matrix, and whose lengths are inversely proportional to 
the eigenvalues of the same matrix. Therefore, a- small eigen- 
value corresponds to a large step component in the direction of 
the associated eigenvector. 

The NLSCIDNT step-cutting procedure works as follows: 

(1) set k = 0; 

(2) if • the current cost is less than the previous Newton- 
Raphson iteration's cost, proceed to the next iteration; 
otherwise, 

(j) set j = k + integer [1 + m/5], where m is the number 
of identified parameters, then set k = j; 

(4) reduce the length of the' components of the steo vector 

in the directions associated with the k smallest eigen- 
values and recompute the step vector. ° 

(o) add this new step vector to the previous Newton-Raphson 
iteration's parameter vector and obtain a new parameter 
vector, compute its associated cost, and go back to 
step 2 . 
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If step cuts have occurred previously (k > 0) and a Newton- 
Raphson iteration is completed successfully without resorting to 
more step cuts, then k is reduced by one. When the lowest level 
of diagnostic printout is flagged, a message is printed to inform 
the user how many eigenvalues are permitted their full size in 
computing the step vector (see Figure 5.1g). 

6.2.3 Primary and Secondary Parameters 

The computer program can identify any or all unknown para- 
meters if enough information is available about these parameters 
from the data. To' identify a large number of parameters, a very 
careful procedure is often required to ensure convergence. In 
the procedure which has been found to be most successful, the 
parameters are divided into two or more groups in decreasing order 
of the effect they have on the rotorcraft response. Initially, only 
the most important parameters are identified, leaving the remaining 
ones fixed at a priori values. Once a reasonable convergence is 
achieved on these parameters, the second group may be added to the 
list of identified parameters. The identification is carried 
out using this new set of parameters, which are identified until 
a reasonable convergence is reached. This procedure is repeated 
until all the parameters are included in the identified set. 

6.2.4 Initial Parameter Estimates 


Good a priori values (e.g., from other wind tunnel data or 
flight tests) should be used for start-up. If they are not avail- 
able, it may be necessary to use a least-squares type of proced- 
ure to obtain starting values. If the start-up values are too far 
from the true answer, the program may converge to a relative mini- 
mum away from the global minimum of the cost function. 
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6.2.5 Data Record Length 


Sufficient data length should be used to identify parameters. 
If the data length is too short, the identified model may yield 
a good time history match even though the parameter estimates 
are inaccurate. A data length equal to 2 to 3 times the longest 
modal period of the system should prove adequate. 

6.2.6 Data Sampling Interval 

In order to assure adequate information content in the data, 
a sampling interval of at least 25 times the highest system 
natural frequency is required. A faster sampling rate provides 
somewhat more accurate parameter estimates because of some addi- 
tional information available in the sampled data. However, the 
algorithm realizes diminishing returns in terms of increased 
parameter accuracy as the sampling rate approaches the continuous 
time case. 

6.3 A PRIORI INFORMATION 

Suppose a user has three data records at the same flight 
condition. Because the control input time histories differ 
among the three records, the information contents of the three 
records are different. By processing these records individually, 
the user will get three slightly different sets of parameter 
estimates and three sets of standard deviations. But what he 
wants is a single set which combines the information content of 
all three records. 

When he has processed the first record and is about to pro-, 
cess the second, he has a priori information from the first 
record about the parameters he is about to identify in the 
second record. He can combine the information of the two records 
by reading into the second NLSCIDNT run: (a) the converged para- 
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meter estimates of the first run as initial parameter estimates 
in the second, and (b) the converged information matrix of the 
first run as the a priori information matrix of the second. 

After the second run has converged, it uses this a priori infor- 
mation with the second record's information to produce a statis- 
tically combined final information matrix and parameter estimates 
according to Eqs . (5.6) and (5.7). 

Similarly, when he processes the third record, he should 
read the combined parameter estimates and information matrix into 
the third run as the initial parameter estimates and a priori 
information matrix. After this run converges, the statistical 
combining is again performed, and the user obtains his goal--one 
set of parameter estimates and standard deviations combining the 
information of all three data records. 

The user must take care that each run identifies exactly 
the same parameters in exactly the same order. The program has 
no way of checking whether this is true. If it is not true, the 
results are invalid. 

The user may find it more convenient to perform this statis- 
tical combining of information himself after processing each of 
the three records individually. In this case, he must write a 
program to do the following: 


wh ere 


n 

M = I M. 
i=l 1 


(6.3) 


P* = M“ a I M.p . 

i-1 1 1 


(6.4) 


Mji is the converged information matrix from the i-th 
run, 

p, is the converged parameter estimates from the i-th 
run , 
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The M 
5 . 2 ). 


n is the number of runs whose answers are combined, 

M is the final information matrix, and 
p* is the combined parameter estimates, 
and p^ may be stored on tape by each run (see Section 
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CHAPTER VII 


DIAGNOSTICS 


7.1 NORMAL TERMINATION PROBLEMS 

Even though a run terminates normally by announcing it has 
converged, there may be problems with the result. Three prob- 
lems are discussed below. 

7.1.1 Large Number of Step-Cuts 

If there were several step-cuts taken just before the search 
"converged,” it may not actually be near a minimum, as explained 
in Section 6 .2 . 2 . Corrective action: 

Cl) Delete parameters having low identifiabilitv from the ■ 
search. 

(2) Restart the . program using the parameter values, from 
the final iteration of the terminated run. 

Increasing the step 7 cut limit (maximum numbe.r_pf. inner iterations per_; 
mitted) beyond 4 or 5 is generally not helpful but may be tried. 

7.1.2 High Value for a Measurement Noise's Variance 

The diagonal elements of the final RCMP matrix are estimates 
of the measurement noises' variances. If one or more are un- 
reasonably large, a problem usually exists. Corrective action: 

(1) Inspect the time history data for errors, particularly 
for disagreements with the program's sign conventions, 
for time lags (for air data errors, particularly), and 
for lost data due to a telemetry dropout. 

( 2 ) Inspect the initial parameter values, particularly for 
incorrect signs or units. 
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(3) Verify that the model equations are reasonable given the 
response data. 2 * * 5 

7.1-3 High Standard Deviation for a Parameter 


If the standard deviation for a particular parameter is high, 
it implies that its identifiability is low. There are several 
possible causes: 

(1) the parameter inherently has little effect on the 
aircraft's response, 

there insufficient excitation of the mode for which 
that parameter is most influential, 

(3) too many parameters are being identified relative to 
the information content of the data record, or 

(4) the parameter's influence is masked by the parallel 
influence of a more important parameter. 

Corrective action: 

( ^ ) -fix the parameter and identify more important para* 
meters, then fix them at their final values and trv 
again to identify the troublesome parameter. 

(2) Try to find a data record having more information on 

the parameter and process it. 


7.2 ABNORMAL TERMINATION PROBLEMS 

Ij. possible, an error message is printed when the program 
terminates abnormally. The conditions resulting in the error 
messages and recommended corrective actions (if not obvious) are 
discussed below. 
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*** ERROR IN SUBROUTINE INPUT. P( j) .NE. 0 BUT AN 

INCORRECT COEFFICIENT INDEX WAS SPECIFIED. 

Parameters 1 through 400 are set aside for use in the poly- 
nomial expansions defining the aerodynamic coefficients. The 
j-th parameter was given a value but not assigned to the expan- 
sion of any aero coefficient ; i.e. , the input card describing the 
j-th parameter (card type 5 in Table 4.1) did not correctly speci 
fy the associated coefficient index number. Valid coefficient 
index numbers are in the range 1 to 384. 

*** ERROR IN SUBROUTINE INPUT. A N0N-ZER0 EXPONENT WAS 

DETECTED F0R AN UNDEFINED EXPANSION VARIABLE. 

If only n(£ 5) expansion variables are defined for the run, 
then values for exponents should not appear for expansion vari- 
ables ^n + l , ‘‘" , ^5 P arameter input cards (type 5 in Table 
4.1). However, one was detected. 

*** ERROR IN SUBROUTINE INPUT. LTYPE = xx, WHICH IS N0T 0NE 

OF THE ALLOWABLE SET: 

Invalid characters were detected in the first two columns 
of a card of type 4 (Table 4.1) . 

*** ERROR IN SUBROUTINE FLAGIT. LL ( j ) = xx, WHICH IS N0T 

IN THE ALLOWABLE SET: 

An ivalid. input name of a state, measurement, or expansion 
variable was detected in an input card of type 4 (Table 4.1). 
Valid input names are listed in Table 4.2.' 
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*** ERR0R IN SUBROUTINE FLAGIT. THE NUMBER OF ELEMENTS 
SPECIFIED BY INPUT IS j, WHICH EXCEEDS THE MAXIMUM ALL0WED , 

k . 

The number of variables specified in a card of type 4 (Table 
4.1) exceeds the maximum number for which the program is dimen- 
sioned. 

*** ERR0R IN SUBROUTINE INREAD. DEVICE ERR0R DETECTED WHILE 
ATTEMPTING T0 READ DATA P0INT k. THE TIME AT THE PREVI0US 
POINT, T(K-l), WAS x. 

Device errors usually result from selecting the wrong density 
or parity for a tape or attempting to read a damaged portion of 
the tape. 

*** ERROR IN SUBROUTINE INREAD. END 0F FILE ENCOUNTERED 
WHILE ATTEMPTING T0 READ DATA POINT k. THE TIME AT. THE 
PREVIOUS POINT, T(K-l), WAS x. 

* 

The number of data sample points specified in input card 
type 2 (Table 4.1) exceeds the number of logical records on the 
time history tape. 

*** ERROR IN SUBROUTINE PRNTP . THE INDEX OF THE COEFFICIENT 
CORRESPONDING TO PARAMETER ( j ) IS ZERO. 

Parameters 1 through 400 must appear in the polynomial expan- 
sion of an aero coefficient. 

*** ERROR IN SUBROUTINE UPDATE. INVALID INPUT WAS PASSED 
TO SUBROUTINE INTGR8. 

The two most likely sources of invalid input to INTGR8 are: 

(1) the maximum number of data points specified in card 
type 2 (Table 4.1) is < 1, or 



(1) restart the program using the parameter values from 
the final iteration of the terminated run, and 

(23 increase the iteration limit. 

*** MAXIMUM STEP-CUTS LIMIT, j, EXCEEDED. PL0TS , IF ANY, 
WILL SH0W Y-HAT FR0M THE PREVI0US ITERATI0N. 

The parameter set used at the beginning of the last Newton- 
Raphson iteration resulted in a higher cost than the previous 
iteration had. The program began cutting back the step in the 
parameters in search of a lower cost, but the number of step-cuts 
exceeded the limit specified in card type 2 (Table 4.1). Correc- 
tive action: 

(1) Delete parameters with low identifiability (small F- 
values) from the search. 

(2) Restart the program Using the parameter values from 
the iteration of the terminated run. 

Increasing the step-cut limit beyond 4 or 5 is generally not 
helpful but may be tried. 

INV ERR0R DETERMINANT OF A=0 RANK OF A=K 

The program attempted to invert a singular matrix. The two 
most frequent causes are a singular inertia matrix in subroutine 
STATE or a. singular information matrix in subroutine SMAIN . 

Both are the result of errors in the program's .input. The infor- 
mation matrix will be singular if, for the vector of measurement 
estimates y and for some one identified parameter 9-, 3y/39.=0 

for the entire time history. For example, if the user mistakenly 
tries to identify the bias in the roll attitude gyro, b 0 , but 

roll attitude is not one of the measurements, then 3v/3b„=0 

0 
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(2) either the relative or absolute error bound specified 
in card type 2 (Table 4.1) is negative, or both are 
zero. 

*** ERR0R IN SUBROUTINE UPDATE. T00 MANY INCREASES IN THE 

INTEGRATION ERROR TOLERANCES. 

If the program cannot integrate the differential equations 
within the error tolerances specified by the user, it increases 
them, prints a message and continues. If this occurs twice, the 
program terminates. Corrective action: examine input error 

tolerances for the possibility of increase or examine data for 
error. 

*** ERROR IN SUBROUTINE UPDATE. M0RE THAN MAXNUM STEPS 

NEEDED BETWEEN T = x AND T* y. 

The program's integration algorithm has determined that in 
order to meet the error bounds it must take more. steps between 
times x and y than the maximum specified in input card type 
3 (Table 4.1). Corrective action: 

(1) verify that the input parameter values are reasonable* 
and if they are, then 

(2) increase the number of steps permitted* or 

(5) relax the error bounds. 

*** EQUATIONS APPEAR T0 BE STIFF 

There are almost certainly errors in the input to the pro- 
gram. However, the differential equations could conceivably 
appear stiff because of exceptionally noisy control variable time 
histories. 
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*** MAXIMUM NUMBER 0F ITERATI0NS EXCEEDED. 

The Levenberg-Marquardt search has not converged in the 
maximum number of iterations specified in card type 2 (Table 4.1). 
Corrective action: re-examine convergence criterion in input 

data and adjust as necessary. 

7.3 ERRORS NOT CAUSING TERMINATION 

Sometimes errors are detected which are not critical to the 
program's computation, and it prints a message and continues. 

The conditions resulting in these error messages are described 
below. 

*** ERROR IN INPUT CARD DECK WHILE ATTEMPTING T0 READ A 
LABEL. ERR0R 0CURRED 0N THIS CARD: 

The character in the first column of an input card of type 
7 (Table 4.1) was invalid. The program skips this card and con- 
tinues to read the next card, which it also expects to be a card 
of type 7 . 

*** INTEGRATION ERR0R TOLERANCES ARE T00 SMALL T0 CONTINUE. 
THEY HAVE 3EEN GIVEN NEW VALUES AT SAMPLE POINT j. RELERR = x 
ABSERR = y. 

The program has determined that it cannot meet the error 
bounds specified by the user and replaces them with the smallest 
bounds it believes it can meet. 

*** ERROR. A NEGATIVE EIGENVALUE, x, WAS COMPUTED FOR THE 
INFORMATION MATRIX. IT HAS BEEN SET = l.E-27. 
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Because the information matrix is theoretically positive 
semidef inite , none of its eigenvalues should be negative. A 
slightly negative eigenvalue occasionally results from numerical 
error when it should actually be slightly positive. A large 
negative eigenvalue invalidates the run, and a serious error 
exists. 

7.4 CONTROL CARDS 

Figure 7.1 shows a typical set of control cards needed to 
run NLSCIDNT . The job and accounting cards are installation 
dependent and so are not shown in detail. The syntax of the 
other control cards may also be installation dependent; the 
cards shown are for the SCOPE 2.1 operating system on the CDC 
7600 computer. 

Tape 2 is the flight test or simulation data to be used 
as input. Tape 3 is the output generated by NLSCIDNT (see 
Section 5.2); if the output does not need to be saved, remove 
the REQUEST and CATALOG cards for tape 3. 
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JOB CARD 
ACCOUNT CARD 

Other cards required by installation to 
initiate a job. 

ATTACH, TAPE 2 = RSRA3, ID = RSH 

REQUEST, TAPE 3, *PF 

ATTACH, NLSCIDNT, ID = BJB 

NLSCIDNT, PL = 7777 

CATALOG, TAPE 3 = RSRAML, ID = RSH. 

7/8/9 (End of Record) 

INPUT DATA CARDS (see Sections 4.1 & 4.2) 

6/7/8/9 (End of File) 


I Figure 7.1 Control Cards 
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APPENDIX A 


MAXIMUM LIKELIHOOD IDENTIFICATION OF PARAMETERS 
OF A NONLINEAR SYSTEM* 


A. 1 INTRODUCTION 

The maximum likelihood method is one of the most flexible 
techniques in statistics for identification of parameters from 
input-output data. Suppose it is possible to make a set of ob- 
servations on a system, whose model has m unknown parameters 0. 
For any given set of values of the parameters 9 from the feasible 
set 0, we can assign a probability p(Z|9) to each outcome Z. 

the outcome of an actual experiment is z, it is of interest 
to know which sets of values of .0 might have led to these obser- 
vations. This concept is embedded in the likelihood function 
. ^(0 | z) . lhis function is of fundamental importance in estimation 
tneory because of tne likelihood principle of Fisher and others 
[6 - 8] which states that if the system model is correct, all in- 
formation about unknown parameters is contained in the likelihood 
runction. The maximum likelihood method finds a set of parameters 
9 to maximize this likelihood function 


0 = max ^(0/z) (A 1 

0 £0 

In other words , the probability of the outcome z is higher with 
parameters 0 in the model than with any other values of para- 
meters from the feasible set. Usually it is more convenient to 
work with the logarithm of the likelihood function (it is possible 
to do so because the logarithm is a strictly monotonic function) . 


This appendix is a shortened version of Appendix A 
Only portions applicable to the NLSCIDNT program 
provide a ready reference for the NLSCIDNT user. 


are 


in Ref. 
included 


1 . 

to 


no 


The great asset of the maximum likelihood method is that it 
can be used with linear or nonlinear models in the presence of 
process and measurement noise. The maximum likelihood estimates 
are asymptotically unbiased, consistent and efficient. These 
terms are defined below. 


Cl) Bias : 0 is an unbiased 


E (9 ) = 0 


where E stands for the 
terms, this implies that 
a large^number of times, 
mates 9 is 9. 


estimator of 9 if 


(A. 2) 


"expected value." In physical 
if the data were collected 
the mean of the different esti- 


(2) ^^•^cient 0 is an efficient estimator of 9 if it is 
unbiased and if for any other unbiased estimator 3 
of 9 


E(0 - 9) 2 < E(9 - 0) 2 


(A. 5) 


(3) Cons is tent : 9 is a consistent estimator of 9 if the 

accuracy of the estimate improves with increasing amount 
of data and approaches the true. value as the amount of 
data tends to infinity. In other words, for anv 
positive n and s ' 


P C I 3 n - 9 | < e) > 1 - n n > N (A 

This definition is analogous to the definition of con- 
vergence in real analysis. 

A. 2 GENERAL NONLINEAR SYSTEMS 


Consider the general nonlinear aircraft equations of motion 
(without process noise) 


x = f(x,u,9,t) 


0 < t < T 


(A. 5) 
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where 

x(t) is n x 1 state vector 

u(t) is Z x 1 input vector 

0 is m x 1 vector of unknown parameters 

Sets of p measurements y(t k ) are taken at discrete times t k 

y (t k ) = h(x(t k ) ,u(t k ) ,9 ,t k ) + v(t k ) 

k - 1,2,3,. ..,N (A • 6 ) 

v(t k ) is Gaussian random noise with the following properties 
E[v(t k )J = 0 

E[v(tj jv^Ctjj)] = R(8,tj)4j k (A. 7) 

The unknown parameters are supposed to occur in the functions 
f and h and in matrices R and x . In the analysis to fol- 
low, the model and the functional form of f and h is assumed 
known correctly. • 

The set of observations y (t^ ,y (t 2 ) , . . . ,y (t N > constitutes 
the outcome 2 in this case. The likelihood function for 0, 

which has the same form as the probability of the outcome z for 

a certain value of parameters 0, is given by 

SB (9|z) = p(zje) 

■ PCyCtj) , y(t 2 ) , . . . , y(t N ) | 9 ) 

= PCY n |9) (A. 8) 
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where 

Y k = ^Ct x ) , . . . ,y(t k ) } , k = 1,2,. ..,N 
PCY N i 0) = pCy(t N )|Y N _ 1 ,0)‘ P CY N _ 1 |3) 

= pCyCt N )|Y N . 1 ,6) P(y(t N . ; )|Y N -2,0) p(Y N _ 2 i3) 

N 

= n pCy(t i )|Y. ,0) (A. 8) 

i»l 

The log- likelihood function is 

N 

log [^(0|i)] = E log (p (y (t • ) | Y. ,,3)} + constant 

i = 1 1 1_i 

(A. 9) 

To find the probability distribution of y(t^) given Y. ^ and 
9, the mean value and covariance are determined first. 

E(yCt i )'| Y i . 1 ,6) A y(i/i-l) (A. 10) 

The expected value or the mean is the best possible estimate of 
measurements at a point given the measurements until the previous 
point. 

cov(y(t i ) lYj.^,6) = H{[y(t i )-y(i/i-l)] [y (t^ -y (i/i-1) ] T } 

A E { v ( i ) v(i) 7 } 

A B Ci) (A. 11) 

v(i) are the innovations at point i and B(i) is the innova- 
tions covariance. Since 
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y(t t ) - ECyCt._)jY i . 1 ,e) - v (i) 


(A. 12) 


it follows that yCt^) given ,Y i _ 1 and 9 have the same distri- 
bution as v(i). Kailath [9] has shown that as the sampling rate 
is increased, the innovations v(i) tend towards having a Gaussian 
density. Assuming a sufficiently high sampling rate, the distri- 
bution of v (i) and, therefore, y(t i ) given Y i _ 1 and 9 is 
Gauss i an , i . e . , 


exp {- f v(i) T B“ 1 Ci)v(i) } 

(yCO I Y- _ 1 , 9 ) = -75 

1 1 r ( 2 | n ■ % | 1/ 2 


(2tt) ' I B (i) 


(A. 15) 


log p(y(.t i ) I Y i _ 1 ,9) 


j v(i) T B' 1 Ci)u(i) 


- j log|B(i)| + constant 


(A.:14) 


The log- likelihood function of Eq. (A. 9). can be written as 

log C 4? (9 | 2 ) ] * - | > Z^{v T (i)B' 1 (i)v(i) + log I B CD t > CA. IS) 

An estimate. of the unknown parameters is obtained by maximizing 
the likelihood function or the log- likelihood function from the 
feasible set of parameter values. 


9 = max log [ 3? (d | z) ] 
9e0 


(A. 16) 


N 


max [- \ Z {v T (i)B' 1 (i)vCi) + log | B(i) | > 1 
9e0 c 1=1 


(A. 17) 
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The log-likelihood function depends on the innovations and 
their covariance. To optimize the likelihood function, a way must 
be found for determining these quantities. In the absence of 
process noise, this is a simple task. 

The state prediction is done using the equations of motion 
(A. 3). The innovation is 

v(i) - y(t i ) - h(J(t i ),u(t i ),9,t i ) . (A. 18) 

From cqs . (A. 6) and (A.”) it can be seen that the covariance of 
the innovations is 

B(i) = R (6 , t i ) (A. 19) 

A. 5 OPTIMIZATION PROCEDURE 

Many possible numerical procedures can be used for this op- 
timization problem. Modified Newton-Raphson [1,2] or Quasilin- 
earization [4] have been found by experience to give quicker con- 
vergence than most procedures like the conjugate gradient or the 
Davison method. The modified Newton-Raphson is a second order 
gradient procedure requiring computation of first and second 
order partials of the log- likelihood function, which, after sub- 
stituting Eq. (A. 19) into (A. 15), is 

log [^C0|2)] = - y L v T (i)R' 1 vCi) + log | R | (A. 20) 

" i=l 
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N 


- - M l p ' 1 m - j x { ^(OR-'Ci) v T (i)R-l(i) 


R ' 1 Ci)vC i J *| Tr^R" 1 Ci) (A. 21) 


o log f g~re I z) i _ 
38 j"k 


N , « T 


2 I^T e " ^ ' R ’ 1(i) 

i=1 i ae k aej 


3v T (i) n-1^ 3R(i) n-L., 

Tet R Cl) “ 38 ~^ R Ci)v(i) 
K J 


3v(i) „-l,,, 3R(i) 

~ Te ' -" ' R ~ R Ci)vu) 


30, 


+ v T (i)R' 1 (i) - R e ^^ R~ 1 (i) - | e Cl) R~ 1 (i)v(i) 


-T Tr 


R^Ci) ^Fr^CI) i|F 

j k J 


*Tr^ R" 1 Ci3vCi) 

j k ' 


-v T (i) R^U) |jS|U R' 1 (i)v (i) 

J ^ 


1 ^)} 


j,k« 1,2,. .. ,p (A. 22) 
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The last three terms in the equation for second partial of the 
log- likelihood function involve second partials of innovation and 
xts covariance. Those terms are dropped because they are much 
smaller than the other terms near a minimum and they are expensive 
to compute. 

Solving Eq. (A. 21) for unknown parameters in R, when R is 
not a function of time, gives 

* 1 N T 

R = N . Z , v v U) (A. 25) 

The equality in (A. 23) holds only for those elements of R which 
are not known a priori. For instance, even if R is known to 
diagonal, the right hand side matrix will not be diagonal in 
general, but the orf-diagonal terms should be ignored before they 
are equated to R. Using (A. 23) in (A. 20) 

N 

log [ (9 j z) ] = - y 2 v T (i) ft ^v(i) + constant (A. 24) 

L i = l 

The first and second derivatives of the log- likelihood func- 
tion with respect to unknown parameters (which do not appear in 
R) are 



log[J?(8|z)] 


N 

V 

A* 

i = l 


v T (i) 


3vCi) 


S 2 log [ ^ (6 j z) ] 

36 j 30 k 


N i 

Z { 
i = l ' 


tv l U) n- 
36. K 


5v(i) 

3 6j 


+ v T (i) 



3 2 v(i) I 
36 ^ 39 ^ ( 


(A. 25) 


(A. 26) 


The terms in the second derivative are approximated as 


3 2 log C6| 2)1 = . 2 / av T (i) ;-i 3v(in 

a6 j 36 k i=l* 39 k "3 / 


(A. 27) 


One Newton-Raphson iteration evaluates the first and second 
gradients for the parameter values given it and computes new para- 
meter values as follows 


0 = 9 , , + A0 

new old 


(A. 28) 


A0 = - 


,2r 

a J 

ae 


)" (»>' 


(A. 29) 


where 


A9 = the step vector, 


J = -log[<?(9|z)b 


3 J 

jq = 1 x m row vector whose j-th element is 3J/39j , and 


3 2 j 

— j = mxm symmetric ^ matrix whose (j,k)th element is 
39 3 2 J/C30 j ae k ) . 


2 9 

Because 3 J/30” is real and symmetric. 


30 


4)" 


^1 T 

• , aT 

1 = 1 1 


(A. 30) 
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where 


Let g 


is the i-th eigenvalue of 3 2 J/ 3 0 ~ , and 
is the i-th eigenvector. 



for simplicity of notation; then 


£0 


m 

2 

i a l 


1 T 

t: Vi* 



(A. 31) 
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